---
title: "Guide to prior setting"
package: MetaboDynamics
author: "Katja Danielzik (katja.danielzik@uni-due.de)"

output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{3. MetaboDynamics_Guide_to_prior_setting}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
set.seed(123)
```

**MetaboDynamics** allows users to specify custom priors for the underlying Bayesian model when using the `model = "raw_plus_counts"` option. This vignette provides a technical guide on how to configure these priors based on your experimental data.

**Note:** If you cannot define appropriate priors for your data, use `model = "scaled_log"`. This option relies on log-transformed, normalized, and scaled cell counts, where priors are fixed to a standard range.

# Prior Parameters

The `raw_plus_counts` model requires three prior specifications. These are passed as arguments to `fit_dynamics_model()`.

| Parameter | Argument | Description |
| :--- | :--- | :--- | :--- |
| **Mean Abundance** | `prior_mean_abundance` | Vector `c(mean, sd)` for the prior on means of log-transformed metabolite abundances. |
| **SD Abundance** |  `prior_sd_abundance` | Scalar representing the *expected value* of the standard deviation of log-transformed metabolite abundances. |
| **Cell Counts** | `prior_counts` | Scalar representing the *expected value* of cell counts. |

**Storage:** If `MetaboDynamics` is performed on a `SummarizedExperiment` object, set priors are stored in `metadata(data)[["priors"]]`.

# Setting Priors Based on Data

To select appropriate values, inspect the distribution of your data. The following examples are simulated data that are similar to our observed experimental data.

## Prior for mean of log-transformed metabolite abundances (`prior_mean_abundance`)

**MetaboDynamics** assumes a metabolite $m$ specific log-normal distribution for 
metabolite abundances. 

\begin{align*}
\log(abundances_{m})&\sim {\sf normal}(\mu_{m},\sigma_{m})
\end{align*}
The prior for $\mu_{m}$ in MetaboDynamics is itself a normal distribution. In the next
I visualize the observed **metabolite specific log-transformed mean abundances** (black)
and the chosen prior (red).

```{r}
library(ggplot2)
library(dplyr)

# Simulate data resembling experimental observations
mu_mean <- 12; sd_mean <- 5
prior_mean_abundance <- c(mu_mean, sd_mean)
prior_sd_abundance <- 2

n_metabolites <- 300
mu <- rnorm(n_metabolites, mean = prior_mean_abundance[1], sd = prior_mean_abundance[2])
sigma <- rexp(n_metabolites, rate = 1/prior_sd_abundance)

n_observations <- 3
abundances <- data.frame(
  metabolite = rep(seq_len(n_metabolites), each = n_observations),
  abundance = NA
)

for (i in seq_len(n_metabolites)){
  abundances[abundances$metabolite == i, ]$abundance <- rnorm(3, mean = mu[i], sd = sigma[i])
}
abundances$abundance <- exp(abundances$abundance)

abundances <- abundances%>%group_by(metabolite)%>%mutate(mean_log=mean(log(abundance)))

# Visualize: Black = Data, Red = Prior
ggplot(abundances, aes(x = mean_log)) +
  geom_density(fill = "grey80", alpha = 0.5) +
  geom_density(aes(x = rnorm(nrow(abundances), mean = 12, sd = 10)), 
               color = "red", linewidth = 1) +
  xlab("log(mean(abundance))")+
  labs(title = "Distribution of metabolite specific mean abundances",
       subtitle = "Black: Observed Data | Red: Prior N(12, 10)")
```

**Recommendation:** Set `prior_mean_abundance` wider than the observed data to
achieve a *weakly informative* prior. In the above plot I chose as a prior for
the mean of log-transformed abundances N(12,10). This translates to the following
parameter setting in `fit_dynamics_model`.
*Example:* `prior_mean_abundance = c(12, 10)`

## Prior for standard deviation of log-transformed metaoblite abundances (`prior_sd_abundance`)

Additional to the prior on the mean $\mu$ of log-transformed metabolite abundances **MetaboDynamics**
also allows to set a prior on the **metabolite** specific standard deviation $\sigma_{mc}$ of log-transformed metabolite abundances `prior_sd_abundance`.

As the observed standard deviation of *metabolite specific* log-transformed metabolite abundances roughly follow an exponential distribution **MetaboDynamics** implements an exponential distribution as the prior for $\sigma$.

**Note**: The *expected value* or mean of the exponential distribution is 
$\frac{1}{rate}$. **MetaboDynamics** takes as input the *expected value* as we
found it more intuitive. If your 
median observed standard deviation is 2, you should set
'prior_sd_abundance' to 2 or greater.

In the next I visualize the observed metabolite specific standard deviation of log-transformed
abundances (black) and the chosen prior (red).

```{r}
sds <- abundances %>%
  group_by(metabolite) %>%
  summarise(sd_log_abundance = sd(log(abundance)))

# Visualize: Black = Data SDs, Red = Prior
ggplot(sds, aes(x = sd_log_abundance)) +
  geom_density(fill = "grey80", alpha = 0.5) +
  geom_density(aes(x = rexp(nrow(sds), rate = 1/4)), # rexp takes as input the rate instead of expected value
               color = "red", linewidth = 1) +
  xlab("sd(log(abundance)")+
  labs(title = "Distribution of metabolite specific standard deviations",
       subtitle = "Black: Observed SDs | Red: Prior Exp(1/4)")
```

**Recommendation:** Set `prior_sd_abundance` slightly larger than the median observed standard deviation, in this way we change the exponential distribution to be wider and therefore the prior to be *weakly informative*. 
*Example:* `prior_sd_abundance = 4`


## Prior for Cell Counts (`prior_counts`)
Similarly to the prior for the standard deviation of log-transformed metabolite abundances **MetaboDynamics**
uses an exponential distribution for cell counts. The argument `prior_counts` sets the *expected value* of cell counts.

**Note**: The *expected value* or mean of the exponential distribution is $\frac{1}{rate}$. **MetaboDynamics** takes as input the *expected value* as we found it more intuitive. 

```{r}
# Simulate counts
counts <- data.frame(
  observation = 1:20, 
  counts = as.integer(rexp(20, 1e-7))
)

# Visualize: Black = Data, Red = Prior
ggplot(counts, aes(x = counts)) +
  geom_density(fill = "grey80", alpha = 0.5) +
  geom_density(data = data.frame(counts = rexp(1e5, 1/5e7)), # rexp takes as input the rate instead of expected value
               aes(x = counts), color = "red", linewidth = 1) +
  labs(title = "Distribution of cell counts",
       subtitle = "Black: Observed Counts | Red: Prior Exp(1/5e7)")
```

**Recommendation:** Set `prior_counts` to a value slightly larger than the median observed counts, in this way we change the exponential distribution to be wider and therefore the prior to be *weakly informative*. 
*Example:* `prior_counts = 5e7`

# Implementation Example

The following code demonstrates how to pass these priors to `fit_dynamics_model()`.

```{r, eval = FALSE}
# library("MetaboDynamics")
# # Prepare data (ensure at least two time points per condition)
# abundances$condition <- "A"
# abundances <- rbind(abundances, abundances)
# abundances$time <- rep(c(1, 2), each = nrow(abundances) / 2)
# 
# counts$condition <- "A"
# counts <- rbind(counts, counts)
# counts$time <- rep(c(1, 2), each = nrow(counts) / 2)
# 
# # Fit model with custom priors
# fit <- fit_dynamics_model(
#   model = "raw_plus_counts",
#   model_option = "sd_per_condition",
#   data = abundances,
#   counts = counts,
#   scaled_measurement = "abundance",
#   
#   # Prior for mean of metabolite specific abundances (Normal)
#   prior_mean_abundance = c(12, 10), 
#   
#   # Prior for metabolite specific standard deviation (Exponential mean)
#   prior_sd_abundance = 4, 
#   
#   # Prior for cell counts (Exponential mean)
#   prior_counts = 5e7
# )
```