# User Manual: IgGeneUsage

SK

#### Sep 12, 2023

# Contents

    require(IgGeneUsage)
    require(rstan)
    require(knitr)
    require(ggplot2)
    require(ggforce)
    require(ggrepel)
    require(reshape2)
    require(patchwork)

# 1 Introduction

Decoding the properties of immune receptor repertoires (IRRs) is key to understanding how our adaptive immune system responds to challenges, such as viral infection or cancer. One important quantitative property of IRRs is their immunoglobulin (Ig) gene usage, i.e. how often are the differnt Igs that make up the immune receptors used in a given IRR. Furthermore, we may ask: is there differential gene usage (DGU) between IRRs from different biological conditions (e.g. healthy vs tumor).

Both of these questions can be answered quantitatively by are answered by IgGeneUsage.

# 2 Input

The main input of IgGeneUsage is a data.frame that has the following columns:

1. individual_id: name of the repertoire (e.g. Patient-1)
2. condition: name of the condition to which each repertoire belongs (healthy, tumor_A, tumor_B, …)
3. gene_name: gene name (e.g. IGHV1-10 or family TRVB1)
4. gene_usage_count: numeric (count) of usage related in individual x gene x condition specified in columns 1-3
5. [optional] repertoire: character/numeric identifier that tags the different biological replicates if they are available for a specific individual

# 3 Model

IgGeneUsage transforms the input data as follows.

First, given \(R\) repertoires with \(G\) genes each, IgGeneUsage generates a gene usage matrix \(Y^{R \times G}\). Row sums in \(Y\) define the total usage (\(N\)) in each repertoire.

Second, for the analysis of DGU between biological conditions, we use a Bayesian model (\(M\)) for zero-inflated beta-binomial regression. Empirically, we know that Ig gene usage data can be noisy also not exhaustive, i.e. some Ig genes that are systematically rearranged at low probability might not be sampled, and certain Ig genes are not encoded (or dysfunctional) in some individuals. \(M\) can fit over-dispersed and zero-inflated Ig gene usage data.

In the output of IgGeneUsage, we report the mean effect size (es or \(\gamma\)) and its 95% highest density interval (HDI). Genes with \(\gamma \neq 0\) (e.g. if 95% HDI of \(\gamma\) excludes 0) are most likely to experience differential usage. Additionally, we report the probability of differential gene usage (\(\pi\)): \[\begin{align} \pi = 2 \cdot \max\left(\int_{\gamma = -\infty}^{0} p(\gamma)\mathrm{d}\gamma, \int_{\gamma = 0}^{\infty} p(\gamma)\mathrm{d}\gamma\right) - 1 \end{align}\] with \(\pi = 1\) for genes with strong differential usage, and \(\pi = 0\) for genes with negligible differential gene usage. Both metrics are computed based on the posterior distribution of \(\gamma\), and are thus related.

# 4 Case Study A: analyzing IRRs

IgGeneUsage has a couple of built-in Ig gene usage datasets. Some were obtained from studies and others were simulated.

Lets look into the simulated dataset `d_zibb_3`. This dataset was generated by a zero-inflated beta-binomial (ZIBB) model, and IgGeneUsage was designed to fit ZIBB-distributed data.

    data("d_zibb_3", package = "IgGeneUsage")
    knitr::kable(head(d_zibb_3))

individual_id | gene_name | gene_usage_count | condition
--- | --- | --- | ---
I_1 | G_1 | 29 | C_1
I_1 | G_2 | 135 | C_1
I_1 | G_3 | 6 | C_1
I_1 | G_4 | 52 | C_1
I_1 | G_5 | 68 | C_1
I_1 | G_6 | 41 | C_1

We can also visualize `d_zibb_3` with ggplot:

    ggplot(data = d_zibb_3)+
      geom_point(aes(x = gene_name, y = gene_usage_count, col = condition),
                 position = position_dodge(width = .7), shape = 21)+
      theme_bw(base_size = 11)+
      ylab(label = "Gene usage [count]")+
      xlab(label = '')+
      theme(legend.position = "top")+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))

## 4.1 DGU analysis

As main input IgGeneUsage uses a data.frame formatted as e.g. `d_zibb_3`. Other input parameters allow you to configure specific settings of the rstan sampler.

In this example, we analyze `d_zibb_3` with 3 MCMC chains, 1500 iterations each including 500 warm-ups using a single CPU core (Hint: for parallel chain execution set parameter `mcmc_cores` = 3). We report for each model parameter its mean and 95% highest density interval (HDIs).

Important remark: you should run DGU analyses using default IgGeneUsage parameters. If warnings or errors are reported with regard to the MCMC sampling, please consult the Stan manual1 1 1 https://mc-stan.org/misc/warnings.html and adjust the inputs accordingly. If the warnings persist, please submit an issue with a reproducible script at the Bioconductor support site or on Github2 2 2 https://github.com/snaketron/IgGeneUsage/issues.

    M <- DGU(ud = d_zibb_3, # input data
             mcmc_warmup = 300, # how many MCMC warm-ups per chain (default: 500)
             mcmc_steps = 1500, # how many MCMC steps per chain (default: 1,500)
             mcmc_chains = 3, # how many MCMC chain to run (default: 4)
             mcmc_cores = 1, # how many PC cores to use? (e.g. parallel chains)
             hdi_lvl = 0.95, # highest density interval level (de fault: 0.95)
             adapt_delta = 0.8, # MCMC target acceptance rate (default: 0.95)
             max_treedepth = 10) # tree depth evaluated at each step (default: 12)

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000493 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.93 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:    1 / 1500 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:   50 / 1500 [  3%]  (Warmup)
    FALSE Chain 1: Iteration:  100 / 1500 [  6%]  (Warmup)
    FALSE Chain 1: Iteration:  150 / 1500 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration:  200 / 1500 [ 13%]  (Warmup)
    FALSE Chain 1: Iteration:  250 / 1500 [ 16%]  (Warmup)
    FALSE Chain 1: Iteration:  300 / 1500 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration:  301 / 1500 [ 20%]  (Sampling)
    FALSE Chain 1: Iteration:  350 / 1500 [ 23%]  (Sampling)
    FALSE Chain 1: Iteration:  400 / 1500 [ 26%]  (Sampling)
    FALSE Chain 1: Iteration:  450 / 1500 [ 30%]  (Sampling)
    FALSE Chain 1: Iteration:  500 / 1500 [ 33%]  (Sampling)
    FALSE Chain 1: Iteration:  550 / 1500 [ 36%]  (Sampling)
    FALSE Chain 1: Iteration:  600 / 1500 [ 40%]  (Sampling)
    FALSE Chain 1: Iteration:  650 / 1500 [ 43%]  (Sampling)
    FALSE Chain 1: Iteration:  700 / 1500 [ 46%]  (Sampling)
    FALSE Chain 1: Iteration:  750 / 1500 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration:  800 / 1500 [ 53%]  (Sampling)
    FALSE Chain 1: Iteration:  850 / 1500 [ 56%]  (Sampling)
    FALSE Chain 1: Iteration:  900 / 1500 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration:  950 / 1500 [ 63%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1500 [ 66%]  (Sampling)
    FALSE Chain 1: Iteration: 1050 / 1500 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 1100 / 1500 [ 73%]  (Sampling)
    FALSE Chain 1: Iteration: 1150 / 1500 [ 76%]  (Sampling)
    FALSE Chain 1: Iteration: 1200 / 1500 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 1250 / 1500 [ 83%]  (Sampling)
    FALSE Chain 1: Iteration: 1300 / 1500 [ 86%]  (Sampling)
    FALSE Chain 1: Iteration: 1350 / 1500 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 1400 / 1500 [ 93%]  (Sampling)
    FALSE Chain 1: Iteration: 1450 / 1500 [ 96%]  (Sampling)
    FALSE Chain 1: Iteration: 1500 / 1500 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 2.865 seconds (Warm-up)
    FALSE Chain 1:                7.245 seconds (Sampling)
    FALSE Chain 1:                10.11 seconds (Total)
    FALSE Chain 1:
    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 2).
    FALSE Chain 2:
    FALSE Chain 2: Gradient evaluation took 0.000162 seconds
    FALSE Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.62 seconds.
    FALSE Chain 2: Adjust your expectations accordingly!
    FALSE Chain 2:
    FALSE Chain 2:
    FALSE Chain 2: Iteration:    1 / 1500 [  0%]  (Warmup)
    FALSE Chain 2: Iteration:   50 / 1500 [  3%]  (Warmup)
    FALSE Chain 2: Iteration:  100 / 1500 [  6%]  (Warmup)
    FALSE Chain 2: Iteration:  150 / 1500 [ 10%]  (Warmup)
    FALSE Chain 2: Iteration:  200 / 1500 [ 13%]  (Warmup)
    FALSE Chain 2: Iteration:  250 / 1500 [ 16%]  (Warmup)
    FALSE Chain 2: Iteration:  300 / 1500 [ 20%]  (Warmup)
    FALSE Chain 2: Iteration:  301 / 1500 [ 20%]  (Sampling)
    FALSE Chain 2: Iteration:  350 / 1500 [ 23%]  (Sampling)
    FALSE Chain 2: Iteration:  400 / 1500 [ 26%]  (Sampling)
    FALSE Chain 2: Iteration:  450 / 1500 [ 30%]  (Sampling)
    FALSE Chain 2: Iteration:  500 / 1500 [ 33%]  (Sampling)
    FALSE Chain 2: Iteration:  550 / 1500 [ 36%]  (Sampling)
    FALSE Chain 2: Iteration:  600 / 1500 [ 40%]  (Sampling)
    FALSE Chain 2: Iteration:  650 / 1500 [ 43%]  (Sampling)
    FALSE Chain 2: Iteration:  700 / 1500 [ 46%]  (Sampling)
    FALSE Chain 2: Iteration:  750 / 1500 [ 50%]  (Sampling)
    FALSE Chain 2: Iteration:  800 / 1500 [ 53%]  (Sampling)
    FALSE Chain 2: Iteration:  850 / 1500 [ 56%]  (Sampling)
    FALSE Chain 2: Iteration:  900 / 1500 [ 60%]  (Sampling)
    FALSE Chain 2: Iteration:  950 / 1500 [ 63%]  (Sampling)
    FALSE Chain 2: Iteration: 1000 / 1500 [ 66%]  (Sampling)
    FALSE Chain 2: Iteration: 1050 / 1500 [ 70%]  (Sampling)
    FALSE Chain 2: Iteration: 1100 / 1500 [ 73%]  (Sampling)
    FALSE Chain 2: Iteration: 1150 / 1500 [ 76%]  (Sampling)
    FALSE Chain 2: Iteration: 1200 / 1500 [ 80%]  (Sampling)
    FALSE Chain 2: Iteration: 1250 / 1500 [ 83%]  (Sampling)
    FALSE Chain 2: Iteration: 1300 / 1500 [ 86%]  (Sampling)
    FALSE Chain 2: Iteration: 1350 / 1500 [ 90%]  (Sampling)
    FALSE Chain 2: Iteration: 1400 / 1500 [ 93%]  (Sampling)
    FALSE Chain 2: Iteration: 1450 / 1500 [ 96%]  (Sampling)
    FALSE Chain 2: Iteration: 1500 / 1500 [100%]  (Sampling)
    FALSE Chain 2:
    FALSE Chain 2:  Elapsed Time: 2.715 seconds (Warm-up)
    FALSE Chain 2:                7.646 seconds (Sampling)
    FALSE Chain 2:                10.361 seconds (Total)
    FALSE Chain 2:
    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 3).
    FALSE Chain 3:
    FALSE Chain 3: Gradient evaluation took 0.000273 seconds
    FALSE Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.73 seconds.
    FALSE Chain 3: Adjust your expectations accordingly!
    FALSE Chain 3:
    FALSE Chain 3:
    FALSE Chain 3: Iteration:    1 / 1500 [  0%]  (Warmup)
    FALSE Chain 3: Iteration:   50 / 1500 [  3%]  (Warmup)
    FALSE Chain 3: Iteration:  100 / 1500 [  6%]  (Warmup)
    FALSE Chain 3: Iteration:  150 / 1500 [ 10%]  (Warmup)
    FALSE Chain 3: Iteration:  200 / 1500 [ 13%]  (Warmup)
    FALSE Chain 3: Iteration:  250 / 1500 [ 16%]  (Warmup)
    FALSE Chain 3: Iteration:  300 / 1500 [ 20%]  (Warmup)
    FALSE Chain 3: Iteration:  301 / 1500 [ 20%]  (Sampling)
    FALSE Chain 3: Iteration:  350 / 1500 [ 23%]  (Sampling)
    FALSE Chain 3: Iteration:  400 / 1500 [ 26%]  (Sampling)
    FALSE Chain 3: Iteration:  450 / 1500 [ 30%]  (Sampling)
    FALSE Chain 3: Iteration:  500 / 1500 [ 33%]  (Sampling)
    FALSE Chain 3: Iteration:  550 / 1500 [ 36%]  (Sampling)
    FALSE Chain 3: Iteration:  600 / 1500 [ 40%]  (Sampling)
    FALSE Chain 3: Iteration:  650 / 1500 [ 43%]  (Sampling)
    FALSE Chain 3: Iteration:  700 / 1500 [ 46%]  (Sampling)
    FALSE Chain 3: Iteration:  750 / 1500 [ 50%]  (Sampling)
    FALSE Chain 3: Iteration:  800 / 1500 [ 53%]  (Sampling)
    FALSE Chain 3: Iteration:  850 / 1500 [ 56%]  (Sampling)
    FALSE Chain 3: Iteration:  900 / 1500 [ 60%]  (Sampling)
    FALSE Chain 3: Iteration:  950 / 1500 [ 63%]  (Sampling)
    FALSE Chain 3: Iteration: 1000 / 1500 [ 66%]  (Sampling)
    FALSE Chain 3: Iteration: 1050 / 1500 [ 70%]  (Sampling)
    FALSE Chain 3: Iteration: 1100 / 1500 [ 73%]  (Sampling)
    FALSE Chain 3: Iteration: 1150 / 1500 [ 76%]  (Sampling)
    FALSE Chain 3: Iteration: 1200 / 1500 [ 80%]  (Sampling)
    FALSE Chain 3: Iteration: 1250 / 1500 [ 83%]  (Sampling)
    FALSE Chain 3: Iteration: 1300 / 1500 [ 86%]  (Sampling)
    FALSE Chain 3: Iteration: 1350 / 1500 [ 90%]  (Sampling)
    FALSE Chain 3: Iteration: 1400 / 1500 [ 93%]  (Sampling)
    FALSE Chain 3: Iteration: 1450 / 1500 [ 96%]  (Sampling)
    FALSE Chain 3: Iteration: 1500 / 1500 [100%]  (Sampling)
    FALSE Chain 3:
    FALSE Chain 3:  Elapsed Time: 2.889 seconds (Warm-up)
    FALSE Chain 3:                7.266 seconds (Sampling)
    FALSE Chain 3:                10.155 seconds (Total)
    FALSE Chain 3:

## 4.2 Output format

In the output of DGU, we provide the following objects:

- `dgu` and `dgu_prob` (main results of IgGeneUsage): quantitative DGU summary on a log- and probability-scale, respectively.
- `gu`: condition-specific relative gene usage (GU) of each gene
- `theta`: probabilities of gene usage in each sample
- `ppc`: posterior predictive checks data (see section ‘Model checking’)
- `ud`: processed Ig gene usage data
- `fit`: rstan (‘stanfit’) object of the fitted model \(\rightarrow\) used for model checks (see section ‘Model checking’)

    summary(M)

    FALSE          Length Class      Mode
    FALSE dgu       9     data.frame list
    FALSE dgu_prob  9     data.frame list
    FALSE gu        8     data.frame list
    FALSE theta    12     data.frame list
    FALSE ppc       2     -none-     list
    FALSE ud       24     -none-     list
    FALSE fit       1     stanfit    S4

## 4.3 Model checking

- Check your model fit. For this, you can use the object glm.

  - Minimal checklist of successful MCMC sampling3 3 3 https://mc-stan.org/misc/warnings.html:
    - no divergences
    - no excessive warnings from rstan
    - Rhat < 1.05
    - high Neff
  - Minimal checklist for valid model:
    - posterior predictive checks (PPCs): is model consistent with reality, i.e. is there overlap between simulated and observed data?
    - leave-one-out analysis

### 4.3.1 MCMC sampling

- divergences, tree-depth, energy
- none found

    rstan::check_hmc_diagnostics(M$fit)

    FALSE
    FALSE Divergences:

    FALSE
    FALSE Tree depth:

    FALSE
    FALSE Energy:

- rhat < 1.05 and n_eff > 0

    rstan::stan_rhat(object = M$fit)|rstan::stan_ess(object = M$fit)

## 4.4 PPC: posterior predictive checks

### 4.4.1 PPCs: repertoire-specific

The model used by IgGeneUsage is generative, i.e. with the model we can generate usage of each Ig gene in a given repertoire (y-axis). Error bars show 95% HDI of mean posterior prediction. The predictions can be compared with the observed data (x-axis). For points near the diagonal \(\rightarrow\) accurate prediction.

    ggplot(data = M$ppc$ppc_rep)+
      facet_wrap(facets = ~individual_id, ncol = 5)+
      geom_abline(intercept = 0, slope = 1, linetype = "dashed", col = "darkgray")+
      geom_errorbar(aes(x = observed_count, y = ppc_mean_count,
                        ymin = ppc_L_count, ymax = ppc_H_count), col = "darkgray")+
      geom_point(aes(x = observed_count, y = ppc_mean_count), size = 1)+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      xlab(label = "Observed usage [counts]")+
      ylab(label = "PPC usage [counts]")

### 4.4.2 PPCs: overall

Prediction of generalized gene usage within a biological condition is also possible. We show the predictions (y-axis) of the model, and compare them against the observed mean usage (x-axis). If the points are near the diagonal \(\rightarrow\) accurate prediction. Errors are 95% HDIs of the mean.

    ggplot(data = M$ppc$ppc_condition)+
      geom_errorbar(aes(x = gene_name, ymin = ppc_L_prop*100,
                        ymax = ppc_H_prop*100, col = condition),
                    position = position_dodge(width = 0.65), width = 0.1)+
      geom_point(aes(x = gene_name, y = ppc_mean_prop*100,col = condition),
                    position = position_dodge(width = 0.65))+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      xlab(label = "Observed usage [%]")+
      ylab(label = "PPC usage [%]")+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))

## 4.5 Results

Each row of `glm` summarizes the degree of DGU observed for specific Igs. Two metrics are reported:

- `es` (also referred to as \(\gamma\)): effect size on DGU, where `contrast` gives the direction of the effect (e.g. tumor - healthy or healthy - tumor)
- `pmax` (also referred to as \(\pi\)): probability of DGU (parameter \(\pi\) from model \(M\))

For `es` we also have the mean, median standard error (se), standard deviation (sd), L (low bound of 95% HDI), H (high bound of 95% HDI)

    kable(x = head(M$dgu), row.names = FALSE, digits = 2)

es_mean | es_mean_se | es_sd | es_median | es_L | es_H | contrast | gene_name | pmax
--- | --- | --- | --- | --- | --- | --- | --- | ---
0.19 | 0.01 | 0.30 | 0.12 | -0.25 | 0.95 | C_1-vs-C_2 | G_1 | 0.48
-0.01 | 0.00 | 0.20 | -0.01 | -0.43 | 0.41 | C_1-vs-C_2 | G_4 | 0.04
-0.09 | 0.00 | 0.23 | -0.06 | -0.60 | 0.34 | C_1-vs-C_2 | G_3 | 0.29
-0.06 | 0.00 | 0.18 | -0.04 | -0.44 | 0.29 | C_1-vs-C_2 | G_2 | 0.26
0.06 | 0.00 | 0.18 | 0.04 | -0.26 | 0.45 | C_1-vs-C_2 | G_5 | 0.27
-0.03 | 0.00 | 0.19 | -0.02 | -0.43 | 0.37 | C_1-vs-C_2 | G_8 | 0.10

### 4.5.1 DGU: differential gene usage

We know that the values of `\gamma` and `\pi` are related to each other. Lets visualize them for all genes (shown as a point). Names are shown for genes associated with \(\pi \geq 0.95\). Dashed horizontal line represents null-effect (\(\gamma = 0\)).

Notice that the gene with \(\pi \approx 1\) also has an effect size whose 95% HDI (error bar) does not overlap the null-effect. The genes with high degree of differential usage are easy to detect with this figure.

    # format data
    stats <- M$dgu
    stats <- stats[order(abs(stats$es_mean), decreasing = FALSE), ]
    stats$gene_fac <- factor(x = stats$gene_name, levels = unique(stats$gene_name))

    ggplot(data = stats)+
      geom_hline(yintercept = 0, linetype = "dashed", col = "gray")+
      geom_errorbar(aes(x = pmax, y = es_mean, ymin = es_L, ymax = es_H),
                    col = "darkgray")+
      geom_point(aes(x = pmax, y = es_mean, col = contrast))+
      geom_text_repel(data = stats[stats$pmax >= 0.95, ],
                      aes(x = pmax, y = es_mean, label = gene_fac),
                      min.segment.length = 0, size = 2.75)+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      xlab(label = expression(pi))+
      xlim(c(0, 1))+
      ylab(expression(gamma))

### 4.5.2 Promising hits

Lets visualize the observed data of the genes with high probability of differential gene usage (\(\pi \geq 0.95\)). Here we show the gene usage in %.

    promising_genes <- stats$gene_name[stats$pmax >= 0.95]

    ppc_gene <- M$ppc$ppc_condition
    ppc_gene <- ppc_gene[ppc_gene$gene_name %in% promising_genes, ]

    ppc_rep <- M$ppc$ppc_rep
    ppc_rep <- ppc_rep[ppc_rep$gene_name %in% promising_genes, ]

    ggplot()+
      geom_point(data = ppc_rep,
                 aes(x = gene_name, y = observed_prop*100, col = condition),
                 size = 1, fill = "black",
                 position = position_jitterdodge(jitter.width = 0.1,
                                                 jitter.height = 0,
                                                 dodge.width = 0.35))+
      geom_errorbar(data = ppc_gene,
                    aes(x = gene_name, ymin = ppc_L_prop*100,
                        ymax = ppc_H_prop*100, group = condition),
                    position = position_dodge(width = 0.35), width = 0.15)+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))+
      ylab(label = "PPC usage [%]")+
      xlab(label = '')

### 4.5.3 Promising hits [count]

Lets also visualize the predicted gene usage counts in each repertoire.

    ggplot()+
      geom_point(data = ppc_rep,
                 aes(x = gene_name, y = observed_count, col = condition),
                 size = 1, fill = "black",
                 position = position_jitterdodge(jitter.width = 0.1,
                                                 jitter.height = 0,
                                                 dodge.width = 0.5))+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      ylab(label = "PPC usage [count]")+
      xlab(label = '')+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))

## 4.6 GU: gene usage summary

IgGeneUsage also reports the inferred gene usage (GU) probability of individual genes in each condition. For a given gene we report its mean GU (`prob_mean`) and the 95% (for instance) HDI (`prob_L` and `prob_H`).

    ggplot(data = M$gu)+
      geom_errorbar(aes(x = gene_name, y = prob_mean, ymin = prob_L,
                        ymax = prob_H, col = condition),
                    width = 0.1, position = position_dodge(width = 0.4))+
      geom_point(aes(x = gene_name, y = prob_mean, col = condition), size = 1,
                 position = position_dodge(width = 0.4))+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      ylab(label = "GU [probability]")+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))

# 5 Leave-one-out (LOO) analysis

To assert the robustness of the probability of DGU (\(\pi\)) and the effect size (\(\gamma\)), IgGeneUsage has a built-in procedure for fully Bayesian leave-one-out (LOO) analysis.

During each step of LOO, we discard the data of one of the R repertoires, and use the remaining data to analyze for DGU. In each step we record \(\pi\) and \(\gamma\) for all genes, including the mean and 95% HDI of \(\gamma\). We assert quantitatively the robustness of \(\pi\) and \(\gamma\) by evaluating their variability for a specific gene.

This analysis can be computationally demanding.

    L <- LOO(ud = d_zibb_3, # input data
             mcmc_warmup = 500, # how many MCMC warm-ups per chain (default: 500)
             mcmc_steps = 1000, # how many MCMC steps per chain (default: 1,500)
             mcmc_chains = 1, # how many MCMC chain to run (default: 4)
             mcmc_cores = 1, # how many PC cores to use? (e.g. parallel chains)
             hdi_lvl = 0.95, # highest density interval level (de fault: 0.95)
             adapt_delta = 0.8, # MCMC target acceptance rate (default: 0.95)
             max_treedepth = 10) # tree depth evaluated at each step (default: 12)

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000172 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.72 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 3.702 seconds (Warm-up)
    FALSE Chain 1:                2.947 seconds (Sampling)
    FALSE Chain 1:                6.649 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.00016 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.6 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 3.862 seconds (Warm-up)
    FALSE Chain 1:                2.544 seconds (Sampling)
    FALSE Chain 1:                6.406 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000167 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.67 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 3.962 seconds (Warm-up)
    FALSE Chain 1:                2.685 seconds (Sampling)
    FALSE Chain 1:                6.647 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000231 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.31 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 3.727 seconds (Warm-up)
    FALSE Chain 1:                2.874 seconds (Sampling)
    FALSE Chain 1:                6.601 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000198 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.98 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.069 seconds (Warm-up)
    FALSE Chain 1:                2.529 seconds (Sampling)
    FALSE Chain 1:                6.598 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000168 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.68 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.164 seconds (Warm-up)
    FALSE Chain 1:                3.257 seconds (Sampling)
    FALSE Chain 1:                7.421 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000159 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.59 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.363 seconds (Warm-up)
    FALSE Chain 1:                2.968 seconds (Sampling)
    FALSE Chain 1:                7.331 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000169 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.69 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.242 seconds (Warm-up)
    FALSE Chain 1:                2.762 seconds (Sampling)
    FALSE Chain 1:                7.004 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000243 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.43 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 3.551 seconds (Warm-up)
    FALSE Chain 1:                2.417 seconds (Sampling)
    FALSE Chain 1:                5.968 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000244 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.44 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.372 seconds (Warm-up)
    FALSE Chain 1:                2.91 seconds (Sampling)
    FALSE Chain 1:                7.282 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000163 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.63 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 3.886 seconds (Warm-up)
    FALSE Chain 1:                2.934 seconds (Sampling)
    FALSE Chain 1:                6.82 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000273 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.73 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.644 seconds (Warm-up)
    FALSE Chain 1:                2.986 seconds (Sampling)
    FALSE Chain 1:                7.63 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000198 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.98 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.154 seconds (Warm-up)
    FALSE Chain 1:                3.345 seconds (Sampling)
    FALSE Chain 1:                7.499 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000234 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.34 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.352 seconds (Warm-up)
    FALSE Chain 1:                3.047 seconds (Sampling)
    FALSE Chain 1:                7.399 seconds (Total)
    FALSE Chain 1:

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000167 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.67 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:  50 / 1000 [  5%]  (Warmup)
    FALSE Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration: 150 / 1000 [ 15%]  (Warmup)
    FALSE Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration: 250 / 1000 [ 25%]  (Warmup)
    FALSE Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration: 350 / 1000 [ 35%]  (Warmup)
    FALSE Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
    FALSE Chain 1: Iteration: 450 / 1000 [ 45%]  (Warmup)
    FALSE Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
    FALSE Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration: 550 / 1000 [ 55%]  (Sampling)
    FALSE Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration: 650 / 1000 [ 65%]  (Sampling)
    FALSE Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 750 / 1000 [ 75%]  (Sampling)
    FALSE Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 850 / 1000 [ 85%]  (Sampling)
    FALSE Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 950 / 1000 [ 95%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 4.565 seconds (Warm-up)
    FALSE Chain 1:                3.062 seconds (Sampling)
    FALSE Chain 1:                7.627 seconds (Total)
    FALSE Chain 1:

Next, we collected the results (GU and DGU) from each LOO iteration:

    L_gu <- do.call(rbind, lapply(X = L, FUN = function(x){return(x$gu)}))
    L_dgu <- do.call(rbind, lapply(X = L, FUN = function(x){return(x$dgu)}))

… and plot them:

## 5.1 LOO-DGU: variability of effect size \(\gamma\)

    ggplot(data = L_dgu)+
      facet_wrap(facets = ~contrast, ncol = 1)+
      geom_hline(yintercept = 0, linetype = "dashed", col = "gray")+
      geom_errorbar(aes(x = gene_name, y = es_mean, ymin = es_L,
                        ymax = es_H, col = contrast, group = loo_id),
                    width = 0.1, position = position_dodge(width = 0.75))+
      geom_point(aes(x = gene_name, y = es_mean, col = contrast,
                     group = loo_id), size = 1,
                 position = position_dodge(width = 0.75))+
      theme_bw(base_size = 11)+
      theme(legend.position = "none")+
      ylab(expression(gamma))

## 5.2 LOO-DGU: variability of \(\pi\)

    ggplot(data = L_dgu)+
      facet_wrap(facets = ~contrast, ncol = 1)+
      geom_point(aes(x = gene_name, y = pmax, col = contrast,
                     group = loo_id), size = 1,
                 position = position_dodge(width = 0.5))+
      theme_bw(base_size = 11)+
      theme(legend.position = "none")+
      ylab(expression(pi))

## 5.3 LOO-GU: variability of the gene usage

    ggplot(data = L_gu)+
      geom_hline(yintercept = 0, linetype = "dashed", col = "gray")+
      geom_errorbar(aes(x = gene_name, y = prob_mean, ymin = prob_L,
                        ymax = prob_H, col = condition,
                        group = interaction(loo_id, condition)),
                    width = 0.1, position = position_dodge(width = 1))+
      geom_point(aes(x = gene_name, y = prob_mean, col = condition,
                     group = interaction(loo_id, condition)), size = 1,
                 position = position_dodge(width = 1))+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      ylab("GU [probability]")+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))

# 6 Case Study B: analyzing IRRs containing biological replicates

    data("d_zibb_4", package = "IgGeneUsage")
    knitr::kable(head(d_zibb_4))

individual_id | condition | gene_name | replicate | gene_usage_count
--- | --- | --- | --- | ---
I_1 | C_1 | G_1 | R_1 | 29
I_1 | C_1 | G_2 | R_1 | 66
I_1 | C_1 | G_3 | R_1 | 285
I_1 | C_1 | G_4 | R_1 | 20
I_1 | C_1 | G_5 | R_1 | 38
I_1 | C_1 | G_6 | R_1 | 709

We can also visualize `d_zibb_4` with ggplot:

    ggplot(data = d_zibb_4)+
      geom_point(aes(x = gene_name, y = gene_usage_count, col = condition,
                     shape = replicate), position = position_dodge(width = 0.8))+
      theme_bw(base_size = 11)+
      ylab(label = "Gene usage [count]")+
      xlab(label = '')+
      theme(legend.position = "top")+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))

## 6.1 Modeling

    M <- DGU(ud = d_zibb_4, # input data
             mcmc_warmup = 500, # how many MCMC warm-ups per chain (default: 500)
             mcmc_steps = 1500, # how many MCMC steps per chain (default: 1,500)
             mcmc_chains = 2, # how many MCMC chain to run (default: 4)
             mcmc_cores = 1, # how many PC cores to use? (e.g. parallel chains)
             hdi_lvl = 0.95, # highest density interval level (de fault: 0.95)
             adapt_delta = 0.8, # MCMC target acceptance rate (default: 0.95)
             max_treedepth = 10) # tree depth evaluated at each step (default: 12)

    FALSE
    FALSE SAMPLING FOR MODEL 'dgu_rep' NOW (CHAIN 1).
    FALSE Chain 1:
    FALSE Chain 1: Gradient evaluation took 0.000606 seconds
    FALSE Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.06 seconds.
    FALSE Chain 1: Adjust your expectations accordingly!
    FALSE Chain 1:
    FALSE Chain 1:
    FALSE Chain 1: Iteration:    1 / 1500 [  0%]  (Warmup)
    FALSE Chain 1: Iteration:   50 / 1500 [  3%]  (Warmup)
    FALSE Chain 1: Iteration:  100 / 1500 [  6%]  (Warmup)
    FALSE Chain 1: Iteration:  150 / 1500 [ 10%]  (Warmup)
    FALSE Chain 1: Iteration:  200 / 1500 [ 13%]  (Warmup)
    FALSE Chain 1: Iteration:  250 / 1500 [ 16%]  (Warmup)
    FALSE Chain 1: Iteration:  300 / 1500 [ 20%]  (Warmup)
    FALSE Chain 1: Iteration:  350 / 1500 [ 23%]  (Warmup)
    FALSE Chain 1: Iteration:  400 / 1500 [ 26%]  (Warmup)
    FALSE Chain 1: Iteration:  450 / 1500 [ 30%]  (Warmup)
    FALSE Chain 1: Iteration:  500 / 1500 [ 33%]  (Warmup)
    FALSE Chain 1: Iteration:  501 / 1500 [ 33%]  (Sampling)
    FALSE Chain 1: Iteration:  550 / 1500 [ 36%]  (Sampling)
    FALSE Chain 1: Iteration:  600 / 1500 [ 40%]  (Sampling)
    FALSE Chain 1: Iteration:  650 / 1500 [ 43%]  (Sampling)
    FALSE Chain 1: Iteration:  700 / 1500 [ 46%]  (Sampling)
    FALSE Chain 1: Iteration:  750 / 1500 [ 50%]  (Sampling)
    FALSE Chain 1: Iteration:  800 / 1500 [ 53%]  (Sampling)
    FALSE Chain 1: Iteration:  850 / 1500 [ 56%]  (Sampling)
    FALSE Chain 1: Iteration:  900 / 1500 [ 60%]  (Sampling)
    FALSE Chain 1: Iteration:  950 / 1500 [ 63%]  (Sampling)
    FALSE Chain 1: Iteration: 1000 / 1500 [ 66%]  (Sampling)
    FALSE Chain 1: Iteration: 1050 / 1500 [ 70%]  (Sampling)
    FALSE Chain 1: Iteration: 1100 / 1500 [ 73%]  (Sampling)
    FALSE Chain 1: Iteration: 1150 / 1500 [ 76%]  (Sampling)
    FALSE Chain 1: Iteration: 1200 / 1500 [ 80%]  (Sampling)
    FALSE Chain 1: Iteration: 1250 / 1500 [ 83%]  (Sampling)
    FALSE Chain 1: Iteration: 1300 / 1500 [ 86%]  (Sampling)
    FALSE Chain 1: Iteration: 1350 / 1500 [ 90%]  (Sampling)
    FALSE Chain 1: Iteration: 1400 / 1500 [ 93%]  (Sampling)
    FALSE Chain 1: Iteration: 1450 / 1500 [ 96%]  (Sampling)
    FALSE Chain 1: Iteration: 1500 / 1500 [100%]  (Sampling)
    FALSE Chain 1:
    FALSE Chain 1:  Elapsed Time: 55.113 seconds (Warm-up)
    FALSE Chain 1:                130.468 seconds (Sampling)
    FALSE Chain 1:                185.581 seconds (Total)
    FALSE Chain 1:
    FALSE
    FALSE SAMPLING FOR MODEL 'dgu_rep' NOW (CHAIN 2).
    FALSE Chain 2:
    FALSE Chain 2: Gradient evaluation took 0.000803 seconds
    FALSE Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 8.03 seconds.
    FALSE Chain 2: Adjust your expectations accordingly!
    FALSE Chain 2:
    FALSE Chain 2:
    FALSE Chain 2: Iteration:    1 / 1500 [  0%]  (Warmup)
    FALSE Chain 2: Iteration:   50 / 1500 [  3%]  (Warmup)
    FALSE Chain 2: Iteration:  100 / 1500 [  6%]  (Warmup)
    FALSE Chain 2: Iteration:  150 / 1500 [ 10%]  (Warmup)
    FALSE Chain 2: Iteration:  200 / 1500 [ 13%]  (Warmup)
    FALSE Chain 2: Iteration:  250 / 1500 [ 16%]  (Warmup)
    FALSE Chain 2: Iteration:  300 / 1500 [ 20%]  (Warmup)
    FALSE Chain 2: Iteration:  350 / 1500 [ 23%]  (Warmup)
    FALSE Chain 2: Iteration:  400 / 1500 [ 26%]  (Warmup)
    FALSE Chain 2: Iteration:  450 / 1500 [ 30%]  (Warmup)
    FALSE Chain 2: Iteration:  500 / 1500 [ 33%]  (Warmup)
    FALSE Chain 2: Iteration:  501 / 1500 [ 33%]  (Sampling)
    FALSE Chain 2: Iteration:  550 / 1500 [ 36%]  (Sampling)
    FALSE Chain 2: Iteration:  600 / 1500 [ 40%]  (Sampling)
    FALSE Chain 2: Iteration:  650 / 1500 [ 43%]  (Sampling)
    FALSE Chain 2: Iteration:  700 / 1500 [ 46%]  (Sampling)
    FALSE Chain 2: Iteration:  750 / 1500 [ 50%]  (Sampling)
    FALSE Chain 2: Iteration:  800 / 1500 [ 53%]  (Sampling)
    FALSE Chain 2: Iteration:  850 / 1500 [ 56%]  (Sampling)
    FALSE Chain 2: Iteration:  900 / 1500 [ 60%]  (Sampling)
    FALSE Chain 2: Iteration:  950 / 1500 [ 63%]  (Sampling)
    FALSE Chain 2: Iteration: 1000 / 1500 [ 66%]  (Sampling)
    FALSE Chain 2: Iteration: 1050 / 1500 [ 70%]  (Sampling)
    FALSE Chain 2: Iteration: 1100 / 1500 [ 73%]  (Sampling)
    FALSE Chain 2: Iteration: 1150 / 1500 [ 76%]  (Sampling)
    FALSE Chain 2: Iteration: 1200 / 1500 [ 80%]  (Sampling)
    FALSE Chain 2: Iteration: 1250 / 1500 [ 83%]  (Sampling)
    FALSE Chain 2: Iteration: 1300 / 1500 [ 86%]  (Sampling)
    FALSE Chain 2: Iteration: 1350 / 1500 [ 90%]  (Sampling)
    FALSE Chain 2: Iteration: 1400 / 1500 [ 93%]  (Sampling)
    FALSE Chain 2: Iteration: 1450 / 1500 [ 96%]  (Sampling)
    FALSE Chain 2: Iteration: 1500 / 1500 [100%]  (Sampling)
    FALSE Chain 2:
    FALSE Chain 2:  Elapsed Time: 62.029 seconds (Warm-up)
    FALSE Chain 2:                126.099 seconds (Sampling)
    FALSE Chain 2:                188.128 seconds (Total)
    FALSE Chain 2:

## 6.2 Posterior predictive checks

    ggplot(data = M$ppc$ppc_rep)+
      facet_wrap(facets = ~individual_id, ncol = 3)+
      geom_abline(intercept = 0, slope = 1, linetype = "dashed", col = "darkgray")+
      geom_errorbar(aes(x = observed_count, y = ppc_mean_count,
                        ymin = ppc_L_count, ymax = ppc_H_count), col = "darkgray")+
      geom_point(aes(x = observed_count, y = ppc_mean_count), size = 1)+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      xlab(label = "Observed usage [counts]")+
      ylab(label = "PPC usage [counts]")

## 6.3 Analysis of estimated effect sizes

The top panel shows the average gene usage (GU) in different biological conditions. The bottom panels shows the differential gene usage (DGU) between pairs of biological conditions.

    g1 <- ggplot(data = M$gu)+
      geom_errorbar(aes(x = gene_name, y = prob_mean, ymin = prob_L,
                        ymax = prob_H, col = condition),
                    width = 0.1, position = position_dodge(width = 0.4))+
      geom_point(aes(x = gene_name, y = prob_mean, col = condition), size = 1,
                 position = position_dodge(width = 0.4))+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      ylab(label = "GU [probability]")+
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4))

    stats <- M$dgu
    stats <- stats[order(abs(stats$es_mean), decreasing = FALSE), ]
    stats$gene_fac <- factor(x = stats$gene_name, levels = unique(stats$gene_name))

    g2 <- ggplot(data = stats)+
      facet_wrap(facets = ~contrast)+
      geom_hline(yintercept = 0, linetype = "dashed", col = "gray")+
      geom_errorbar(aes(x = pmax, y = es_mean, ymin = es_L, ymax = es_H),
                    col = "darkgray")+
      geom_point(aes(x = pmax, y = es_mean, col = contrast))+
      geom_text_repel(data = stats[stats$pmax >= 0.95, ],
                      aes(x = pmax, y = es_mean, label = gene_fac),
                      min.segment.length = 0, size = 2.75)+
      theme_bw(base_size = 11)+
      theme(legend.position = "top")+
      xlab(label = expression(pi))+
      xlim(c(0, 1))+
      ylab(expression(gamma))

    (g1/g2)

# 7 Session

    sessionInfo()

    FALSE R version 4.5.0 (2025-04-11 ucrt)
    FALSE Platform: x86_64-w64-mingw32/x64
    FALSE Running under: Windows Server 2022 x64 (build 20348)
    FALSE
    FALSE Matrix products: default
    FALSE   LAPACK version 3.12.1
    FALSE
    FALSE locale:
    FALSE [1] LC_COLLATE=C
    FALSE [2] LC_CTYPE=English_United States.utf8
    FALSE [3] LC_MONETARY=English_United States.utf8
    FALSE [4] LC_NUMERIC=C
    FALSE [5] LC_TIME=English_United States.utf8
    FALSE
    FALSE time zone: America/New_York
    FALSE tzcode source: internal
    FALSE
    FALSE attached base packages:
    FALSE [1] stats     graphics  grDevices utils     datasets  methods   base
    FALSE
    FALSE other attached packages:
    FALSE  [1] patchwork_1.3.0     reshape2_1.4.4      ggrepel_0.9.6
    FALSE  [4] ggforce_0.4.2       ggplot2_3.5.2       knitr_1.50
    FALSE  [7] rstan_2.32.7        StanHeaders_2.32.10 IgGeneUsage_1.23.0
    FALSE [10] BiocStyle_2.37.0
    FALSE
    FALSE loaded via a namespace (and not attached):
    FALSE  [1] SummarizedExperiment_1.39.0 gtable_0.3.6
    FALSE  [3] xfun_0.52                   bslib_0.9.0
    FALSE  [5] QuickJSR_1.7.0              inline_0.3.21
    FALSE  [7] Biobase_2.69.0              lattice_0.22-7
    FALSE  [9] vctrs_0.6.5                 tools_4.5.0
    FALSE [11] generics_0.1.4              stats4_4.5.0
    FALSE [13] curl_6.2.3                  parallel_4.5.0
    FALSE [15] tibble_3.2.1                pkgconfig_2.0.3
    FALSE [17] Matrix_1.7-3                RColorBrewer_1.1-3
    FALSE [19] S4Vectors_0.47.0            RcppParallel_5.1.10
    FALSE [21] lifecycle_1.0.4             stringr_1.5.1
    FALSE [23] compiler_4.5.0              farver_2.1.2
    FALSE [25] tinytex_0.57                codetools_0.2-20
    FALSE [27] GenomeInfoDb_1.45.4         htmltools_0.5.8.1
    FALSE [29] sass_0.4.10                 yaml_2.3.10
    FALSE [31] tidyr_1.3.1                 pillar_1.10.2
    FALSE [33] crayon_1.5.3                jquerylib_0.1.4
    FALSE [35] MASS_7.3-65                 DelayedArray_0.35.1
    FALSE [37] cachem_1.1.0                magick_2.8.6
    FALSE [39] abind_1.4-8                 tidyselect_1.2.1
    FALSE [41] digest_0.6.37               stringi_1.8.7
    FALSE [43] purrr_1.0.4                 dplyr_1.1.4
    FALSE [45] bookdown_0.43               labeling_0.4.3
    FALSE [47] polyclip_1.10-7             fastmap_1.2.0
    FALSE [49] grid_4.5.0                  colorspace_2.1-1
    FALSE [51] cli_3.6.5                   SparseArray_1.9.0
    FALSE [53] magrittr_2.0.3              loo_2.8.0
    FALSE [55] S4Arrays_1.9.1              dichromat_2.0-0.1
    FALSE [57] pkgbuild_1.4.8              withr_3.0.2
    FALSE [59] scales_1.4.0                UCSC.utils_1.5.0
    FALSE [61] rmarkdown_2.29              XVector_0.49.0
    FALSE [63] httr_1.4.7                  matrixStats_1.5.0
    FALSE [65] gridExtra_2.3               evaluate_1.0.3
    FALSE [67] GenomicRanges_1.61.0        IRanges_2.43.0
    FALSE [69] V8_6.0.4                    rstantools_2.4.0
    FALSE [71] rlang_1.1.6                 Rcpp_1.0.14
    FALSE [73] glue_1.8.0                  tweenr_2.0.3
    FALSE [75] BiocManager_1.30.25         BiocGenerics_0.55.0
    FALSE [77] jsonlite_2.0.0              plyr_1.8.9
    FALSE [79] R6_2.6.1                    MatrixGenerics_1.21.0