manual_fitting_species.R

library(tidyverse)
library(caret)
library(fitdistrplus)
library(logspline)
library(goftest)
library(rmarkdown)
library(Metrics)

fit_stable_dist <- function(data, dist) {
  result <- NULL
  
  # Try MLE method
  result <- tryCatch({
    fitdist(data, dist, method = "mle")
  }, error = function(e) {
    NULL
  })
  
  # If MLE fails, try MGE method
  if (is.null(result)) {
    result <- tryCatch({
      fitdist(data, dist, method = "mge")
    }, error = function(e) {
      NULL
    })
  }
  
  # If MGE fails, try MME method
  if (is.null(result)) {
    result <- tryCatch({
      fitdist(data, dist, method = "mme")
    }, error = function(e) {
      NULL
    })
  }
  return(result)
}



fractional_subset <- function(vec, frac) {
  num_elements <- ceiling(length(vec) * frac)
  return(sample(vec, num_elements))
}

all_cpms <- numeric()

otu_numbers <- "/mnt/tlin/vol/jlab/tlin/in_silico/otu_assignment/vsearch/healthy-fecal-table-dn-97.tsv"

otu_df <- read_tsv(otu_numbers, skip = 1)

otu_df <- dplyr::select(otu_df, -c(`#OTU ID`))

filtered_otus <- otu_df %>%
  rowwise() %>%
  mutate(max_value=max(c_across(everything()))) %>%
  filter(max_value > 5) %>%
  as.matrix() %>%
  as.tibble()

#non 0 values per column 

filtered_otus <- otu_df %>%
  rowwise() %>%
  mutate(max_value=max(c_across(everything()))) %>%
  mutate(occurrence=sum(c_across(everything()))!=0) %>%
  filter(max_value > 5 && occurrence > 3) %>%
  as.matrix() %>%
  as.tibble()

filtered_otus <- dplyr::select(filtered_otus, -c(`max_value`))

all_species_counts <- c()

for (column in names(filtered_otus)) {
  all_species_counts <- c(all_species_counts, filtered_otus[[column]])
}

# Convert to a data frame for fitting
num_reads <- data.frame(SpeciesCount = all_species_counts)

positive_num_reads <- all_species_counts[all_species_counts>0]

rm(all_groups)

all_groups <- data.frame(groups = positive_num_reads)


descdist(all_groups$groups, discrete = FALSE)
fg.mle <-fit_stable_dist(all_groups$groups,"gamma")

process <- preProcess(all_groups, method=c("range"))
normalized <- predict(process, all_groups)
normalized_data <- normalized$groups

fit.poisson <- fit_stable_dist(all_groups$groups, "pois")

# Fit Binomial Distribution (assuming 10 trials)
size <- 30
prob_initial <- mean(all_groups$groups) / size

# Fit Binomial Distribution with initial guess for prob
#fit.binomial <- fitdist(all_groups$groups, "binom", fix.arg = list(size = size), start = list(prob = prob_initial))

fit.geometric <- fit_stable_dist(all_groups$groups, "geom")

# Fit Negative Binomial Distribution
fit.negbinom <- fit_stable_dist(all_groups$groups, "nbinom")

fit.normal <- fit_stable_dist(all_groups$groups, "norm")

# Fit Gamma Distribution
fit.gamma <- fit_stable_dist(all_groups$groups, "gamma")

# Fit Exponential Distribution
fit.exponential <- fit_stable_dist(all_groups$groups, "exp")

# Fit Log-Normal Distribution
fit.lognormal <- fit_stable_dist(all_groups$groups, "lnorm")

fit.beta <- fit_stable_dist(normalized_data, "beta")


# Perform Anderson-Darling test for each distribution
ad_poisson <- ad.test(all_groups$groups, ppois, lambda = fit.poisson$estimate)
ad_binomial <- NULL#ad.test(all_groups$groups, pbinom, size = size, prob = fit.binomial$estimate)
ad_geometric <- ad.test(all_groups$groups, pgeom, prob = fit.geometric$estimate)
ad_negbinom <- ad.test(all_groups$groups, pnbinom, size = fit.negbinom$estimate[1], mu = fit.negbinom$estimate[2])
ad_normal <- ad.test(all_groups$groups, pnorm, mean = fit.normal$estimate[1], sd = fit.normal$estimate[2])
ad_gamma <- ad.test(all_groups$groups, pgamma, shape = fit.gamma$estimate[1], rate = fit.gamma$estimate[2])
ad_exponential <- ad.test(all_groups$groups, pexp, rate = fit.exponential$estimate)
ad_lognormal <- ad.test(all_groups$groups, plnorm, meanlog = fit.lognormal$estimate[1], sdlog = fit.lognormal$estimate[2])
ad_beta <- ad.test(normalized_data, pbeta, shape1 = fit.beta$estimate[1], shape2 = fit.beta$estimate[2])

# Save fitted models and additional statistics into the environment for R Markdown
save(fit.poisson, fit.binomial, fit.geometric, fit.negbinom, fit.normal, fit.gamma, fit.exponential, fit.lognormal, fit.beta,
     ad_poisson, ad_binomial, ad_geometric, ad_negbinom, ad_normal, ad_gamma, ad_exponential, ad_lognormal, ad_beta,
     file = "fit_models_spcies_counts.RData")

# Render the R Markdown file to HTML
render("template_curves_cpms.Rmd", output_file = "fit_summaries_and_plots_species_counts.html")

df <- data.frame(Model = character(), MSE = numeric(), AD_Test = numeric(), AD_PVALUE = numeric(), Estimate = character(), AIC = numeric(), stringsAsFactors = FALSE)


distributions <- c("gamma",
                   "pois",
#                   "binom",
                   "geom",
                   "nbinom",
                   "norm",
                   "exp",
                   "lnorm",
                   "weibull",
                   "beta")

# Loop to add rows to the dataframe
#redo the calculation in the loop use a string list for making it easier
for (dist in distributions) {
  print(dist)
  if (dist == "binom") {
    # Fit Binomial Distribution (assuming 10 trials)
    prob_initial <- mean(all_groups$groups)
    
    fit <- fit_stable_dist(all_groups$groups, "binom")
  } else if (dist == "beta") {
    process <- preProcess(all_groups, method=c("range"))
    normalized <- predict(process, all_groups)
    normalized_data <- normalized$groups
    fit <- fit_stable_dist(normalized_data, "beta")
    #ad_beta <- ad.test(normalized_data, pbeta, shape1 = fit.beta$estimate[1], shape2 = fit.beta$estimate[2])
  } else {
    fit <- fit_stable_dist(all_groups$groups, dist)
  }
  print(fit$estimate)
  data <- all_groups$groups
  ad <- switch(dist,
                 gamma = NULL, #ad.test(data, pgamma, shape = fit$estimate[1], rate = fit$estimate[2]),
                 pois = ad.test(data, ppois, lambda = fit$estimate),
                 binom = ad.test(data, pbinom, size = size, prob = fit$estimate),
                 geom = ad.test(data, pgeom, prob = fit$estimate),
                 nbinom = ad.test(data, pnbinom, size = fit$estimate[1], prob = fit$estimate[2] / (fit$estimate[2] + fit$estimate[1])),
                 norm = ad.test(data, pnorm, mean = fit$estimate[1], sd = fit$estimate[2]),
                 exp = NULL, #ad.test(data, pexp, rate = fit$estimate),
                 lnorm = ad.test(data, plnorm, meanlog = fit$estimate[1], sdlog = fit$estimate[2]),
                 beta = ad.test(normalized_data, pbeta, shape1 = fit$estimate[1], shape2 = fit$estimate[2]),
                 weibull = ad.test(data, pweibull, shape = fit$estimate[1], scale = fit$estimate[2]),
                 stop("Unknown distribution"))
  gof <- gofstat(fit)
  # Create a new row as a dataframe
  new_row <- data.frame(Model = dist, 
                        AIC = fit$aic, stringsAsFactors = FALSE,
                        AD_Test = if (!is.null(ad$statistic)) ad$statistic else NA,
                        AD_PVALUE = if (!is.null(ad$p.value)) ad$p.value else NA,
                        Estimate = paste(names(fit$estimate), fit$estimate, sep = "=", collapse = ", "),
                        KS = if (!is.null(gof$ks)) gof$ks else NA,
                        KS_PVALUE = if (!is.null(gof$kstest)) gof$kstest else NA,
                        CvM = if (!is.null(gof$stat)) gof$stat else NA,
                        BIC = if (!is.null(gof$bic)) gof$bic else NA,
                        ChiSq = if (!is.null(gof$chisq)) gof$chisq else NA
  )
                        
  # Add the new row to the dataframe
  df <- rbind(df, new_row)
}

# Display the dataframe
print(df)
write.csv(df, file = "manual_fitted_species_counts.csv")