david-barnett / microViz

R package for microbiome data visualization and statistics. Uses phyloseq, vegan and the tidyverse. Docker image available.
https://david-barnett.github.io/microViz/
GNU General Public License v3.0
94 stars 10 forks source link

Mixed models using `taxatree_models()` #51

Open fconstancias opened 2 years ago

fconstancias commented 2 years ago

Hi @david-barnett,

Thanks again for your effort in developping the package. How would you recommend to run mixed models thourgh taxatree_models?

Following your tutorial, I am not able specify type = lmer

lm_models <- phylo %>% 
  tax_fix() %>% 
  tax_prepend_ranks() %>% 
  # it makes sense to perform the compositional transformation BEFORE filtering
  tax_transform("compositional", rank = "Genus", keep_counts = TRUE) %>% 
  tax_filter(min_prevalence = 0.1, undetected = 0, use_counts = TRUE) %>% 
  tax_transform(
    trans = "log2", chain = TRUE, zero_replace = "halfmin"
  ) %>% 
  taxatree_models(
    type = lmer, 
    ranks = NULL, # uses every rank available except the first
    # variables = c("UC", "female", "antibiotics", "steroids", "age_scaled")
    formula = ~ UC    + (1| female)

  )
Proportional min_prevalence given: 0.1 --> min 7/67 samples.
2022-06-02 10:33:27 - modelling at rank: Phylum
Error in `[[<-`(`*tmp*`, "call", value = f) : 
  [[<- defined for objects of type "S4" only for subclasses of environment

Using the R random effect syntax and type = lm does not generate any issue.

lm_models <- phylo %>% 
  tax_fix() %>% 
  tax_prepend_ranks() %>% 
  # it makes sense to perform the compositional transformation BEFORE filtering
  tax_transform("compositional", rank = "Genus", keep_counts = TRUE) %>% 
  tax_filter(min_prevalence = 0.1, undetected = 0, use_counts = TRUE) %>% 
  tax_transform(
    trans = "log2", chain = TRUE, zero_replace = "halfmin"
  ) %>% 
  taxatree_models(
    type = lm, 
    ranks = NULL, # uses every rank available except the first
    # variables = c("UC", "female", "antibiotics", "steroids", "age_scaled")
    formula = ~ UC    + (1| female)
  )

However, the random effect is not evaluated - which is not surprising.

lm_models %>% 
  taxatree_models2stats() %>% 
  .["taxatree_stats"]

$taxatree_stats
# A tibble: 268 × 7
   term           taxon            rank  estimate std.error statistic  p.value
   <fct>          <chr>            <fct>    <dbl>     <dbl>     <dbl>    <dbl>
 1 UC             P: Firmicutes    Phyl…   -47.4     10.5       -4.50  2.84e-5
 2 1 | femaleTRUE P: Firmicutes    Phyl…    NA       NA         NA    NA      
 3 UC             P: Bacteroidetes Phyl…   -16.7      2.92      -5.72  2.97e-7
 4 1 | femaleTRUE P: Bacteroidetes Phyl…    NA       NA         NA    NA      

Thanks.

david-barnett commented 2 years ago

Hi Florentin!

I fixed this particular bug in a new branch https://github.com/david-barnett/microViz/tree/tax-model-lmer-fix Feel free to install from there and play around, you'll probably need to load broom.mixed or even explicitly provide a tidier function for the lmer models in the taxatree_models2stats fun argument.

I'm reluctant to push this fix to the main branch at the moment as I have not tested out mixed models with tax_model and taxatree_models and don't want to encourage its general use before I do so.

I likely will do this in the next month or so, but I would be grateful to hear more from you if you get this working or encounter further problems 🙂

fconstancias commented 2 years ago

Hi @david-barnett,

Thanks a lot. I will give it a try when I have a bit of time and will keep you updated.

Best