corncob icon indicating copy to clipboard operation
corncob copied to clipboard

Feature request sent via email: extracting estimates

Open bryandmartin opened this issue 4 years ago • 2 comments

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

bryandmartin avatar Oct 22 '20 18:10 bryandmartin

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


ianartmor avatar Nov 04 '20 23:11 ianartmor

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(.)
    )  
}

roey-angel avatar Mar 01 '21 17:03 roey-angel