Contents

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.

1 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"]].

2 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.

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

MetaboDynamics assumes a metabolite 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\) in MetaboDynamics is a normal distribution. In the next I visualize the observed data (black) and the chosen prior (red).

library(ggplot2)
library(dplyr)
library(MetaboDynamics)

# 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)

# Visualize: Black = Data, Red = Prior
ggplot(abundances, aes(x = log(abundance))) +
  geom_density(fill = "grey80", alpha = 0.5) +
  geom_density(aes(x = rnorm(nrow(abundances), mean = 12, sd = 10)), 
               color = "red", linewidth = 1) +
  labs(title = "Distribution of log-transformed metabolite 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)

2.2 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 standard deviation \(\sigma\) 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\). prior_sd_abundanceis hereby the rate (\(\lambda\)) of the exponential distribution of standard deviations.

Note: The expected value/mean of the exponential distribution is 1/rate. MetaboDynamics takes as input the rate (\(\lambda\)) as we found it more intuitive. If your median observed standard deviation is 2, the rate is 2 and the mean equals 1/2.

In the next I visualize the observed data (black) and the chosen prior (red).

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)), 
               color = "red", linewidth = 1) +
  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 counts, in this way we change the exponential distribution to be wider and therefore the prior to be weakly informative. Example: prior_sd_abundance = 4

2.3 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 1/rate. MetaboDynamics takes as input the rate (\(\lambda\)) as we found it more intuitive.

# 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)), 
               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

3 Implementation Example

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

# # 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
# )