In kylebittinger/polyafit: Fit observed data to a multivariate Polya distribution

library(tidyverse)

as_proportion <- function (x) x / sum(x)
set.seed(1)
example_props <- as_proportion(rchisq(50, 2))

set.seed(5)
polyafit_example_data <- t(cbind(
  rmultinom(1, 500, example_props),
  rmultinom(1, 600, example_props),
  rmultinom(1, 550, example_props)))
polyafit_example_data[3,4] <- 150
rownames(polyafit_example_data) <- c("sample_1", "sample_2", "sample_3")
colnames(polyafit_example_data) <- paste("species", 1:50, sep = "_")

polyafit_example_data

fit12 <- polyafit::pfit(polyafit_example_data[c("sample_1", "sample_2"),])

fit12 %>%
  plot()

fit12 %>%
  polyafit::feature_enrichment() %>%
  filter(p.value < 0.05)

fit12$theta

fit13 <- polyafit::pfit(polyafit_example_data[c("sample_1", "sample_3"),])

fit13 %>%
  plot()

fit13 %>%
  polyafit::feature_enrichment() %>%
  filter(p.value < 0.05)

fit13$theta

usethis::use_data(polyafit_example_data, overwrite = TRUE)

kylebittinger/polyafit documentation built on Jan. 11, 2023, 8:53 a.m.