complex mean-variability association

[stemangiola]

Hello @lel4011

FYI we see cases where the intercept association is positive while the slope is negative

In particular, https://github.com/lel4011/sccomp/blob/b46f806a20022edb1ddefb6ff6e37c7ab77c4120/inst/stan/glm_multi_beta_binomial.stan#L568 the priors should be centred on 0 to avoid bias.

I would be curious if we change those priors, what different fit we might get. As the one now is apparently not great

The code used for analysis

suppressPackageStartupMessages({
  library(curatedMetagenomicData)
  library(TreeSummarizedExperiment)
  library(dplyr)
  library(tidyr)
  library(sccomp)
  library(tidySummarizedExperiment)
  library(magrittr)
})

tse <- curatedMetagenomicData(
  "BritoIL_2016.relative_abundance",
  dryrun = FALSE,
  counts = TRUE
) %$%
 `2021-03-31.BritoIL_2016.relative_abundance` |>
  dplyr::filter(body_site %in% c("stool", "oralcavity")) |>
  dplyr::mutate(body_site = factor(body_site, levels = c("stool", "oralcavity")))

# MIA-style feature filtering:
# - detection threshold: at least 10 counts to call a taxon "present"
# - prevalence threshold: present in at least 10% of samples

n_taxa_before <- nrow(tse)
tse <- mia::subsetByPrevalent(
  tse,
  detection = 10L,
  prevalence = 0.10,
  assay.type = "relative_abundance"
)

fit <- mia::meltAssay(
    tse,
    assay.type = "relative_abundance",
    add.row = FALSE,
    add.col = TRUE
  ) |>

  sccomp::sccomp_estimate(
    formula_composition = ~ body_site,
    formula_variability = ~ body_site,
    sample = "SampleID",
    cell_group = "FeatureID",
    abundance = "relative_abundance",
    inference_method = "hmc",
    cores = max(1, parallelly::availableCores() - 1), 
    bimodal_mean_variability_association = TRUE,
    verbose = FALSE, 
    max_sampling_iterations = 4000,
    sig_figs = 6
  )