statdivlab / corncob

Count Regression for Correlated Observations with the Beta-binomial
103 stars 22 forks source link

Feature request sent via email: extracting estimates #93

Open bryandmartin opened 4 years ago

bryandmartin commented 4 years ago

Extract the estimates from differentialTest without having to use a long pipe with map

ianartmor commented 4 years ago

I have a function for this too, so I'll post it in case someone needs it before this feature is added. Like the interval function I posted yesterday, it does require purrr.

extractCoefficients <- function(diftest_res, taxa){

#add names to all_models list
  names(diftest_res$all_models) <- taxa_names(diftest_res$data)

  # default is all taxa
  if (missing(taxa)) {
    taxa <- taxa_names(diftest_res$data)
  } else {
    taxa <- taxa
  }

  purrr::map(
    diftest_res$all_models[taxa],
    ~ stats::coef(.)
  )

}

Example usage:


data("soil_phylo")

soil_diftest_res <- differentialTest(
  data = soil_phylo %>% tax_glom("Genus"), formula = ~DayAmdmt, phi.formula = ~1,
  formula_null = ~1, phi.formula_null = ~1, test = "LRT", full_output = TRUE
)

diftest_coefs <- extractCoefficients(soil_diftest_res, soil_diftest_res$significant_taxa)

#coefficients for first significant model
diftest_coefs[1] 

#can also subset using a character vector because it's a named list
diftest_coefs["OTU.2"]
roey-angel commented 3 years ago

That's a useful function @ianartmor but note that differentialTest() occasionally produces NA for some OTUs and then calling coef() throws an error.

Here's a modified version that works around it but might produce a list that's shorter than the number of original OTUs

extractCoefficients <- function(diftest_res, taxa){

#add names to all_models list
  names(diftest_res$all_models) <- taxa_names(diftest_res$data)

  # default is all taxa
  if (missing(taxa)) {
    taxa <- taxa_names(diftest_res$data)
  } else {
    taxa <- taxa
  }

  diftest_res$all_models %>% 
    map(~ .x[!is.na(.x)]) %>% 
    compact() %>% 
    purrr::map(~ stats::coef(.)
    )  
}