Skip to content

Add example of one sample #108

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
stemangiola wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Add example of one sample #108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6408a14
improve messaging
stemangiola Dec 14, 2023
c615938
changed arguments and add README
stemangiola Dec 15, 2023
2962584
add figures README
stemangiola Dec 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ Imports:
readr,
scales,
stringr,
glue
glue,
fitdistrplus
Suggests:
BiocStyle,
testthat (>= 3.0.0),
Expand Down
1 change: 1 addition & 0 deletions NAMESPACE
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ importFrom(dplyr,left_join)
importFrom(dplyr,mutate)
importFrom(dplyr,select)
importFrom(dplyr,with_groups)
importFrom(fitdistrplus,fitdist)
importFrom(forcats,fct_relevel)
importFrom(forcats,fct_reorder)
importFrom(ggrepel,geom_text_repel)
Expand Down
122 changes: 102 additions & 20 deletions R/functions_multi_beta_binomial.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,8 @@ sccomp_glm_data_frame_raw = function(.data,

# Secondary arguments
contrasts = NULL,
prior_mean = list(intercept = c(0,1), coefficients = c(0,1)),
prior_overdispersion_mean_association = list(intercept = c(5, 2), slope = c(0, 0.6), standard_deviation = c(20, 40)),
prior_mean = list(intercept_normal = c(0,1), coefficients_normal = c(0,1)),
prior_overdispersion_mean_association = list(intercept_normal = c(5, 2), slope_normal = c(0, 0.6), sd_gamma = c(20, 40)),
percent_false_positive = 5,
check_outliers = TRUE,
approximate_posterior_inference = "none",
Expand Down Expand Up @@ -110,8 +110,8 @@ sccomp_glm_data_frame_counts = function(.data,
# Secondary arguments
contrasts = NULL,
#.grouping_for_random_intercept = NULL,
prior_mean = list(intercept = c(0,1), coefficients = c(0,1)),
prior_overdispersion_mean_association = list(intercept = c(5, 2), slope = c(0, 0.6), standard_deviation = c(20, 40)),
prior_mean = list(intercept_normal = c(0,1), coefficients_normal = c(0,1)),
prior_overdispersion_mean_association = list(intercept_normal = c(5, 2), slope_normal = c(0, 0.6), sd_gamma = c(20, 40)),
percent_false_positive = 5,
check_outliers = TRUE,
approximate_posterior_inference = "none",
Expand Down Expand Up @@ -261,11 +261,11 @@ sccomp_glm_data_frame_counts = function(.data,
data_for_model$TNS = length(data_for_model$truncation_not_idx)

# Prior
data_for_model$prior_prec_intercept = prior_overdispersion_mean_association$intercept
data_for_model$prior_prec_slope = prior_overdispersion_mean_association$slope
data_for_model$prior_prec_sd = prior_overdispersion_mean_association$standard_deviation
data_for_model$prior_mean_intercept = prior_mean$intercept
data_for_model$prior_mean_coefficients = prior_mean$coefficients
data_for_model$prior_prec_intercept = prior_overdispersion_mean_association$intercept_normal
data_for_model$prior_prec_slope = prior_overdispersion_mean_association$slope_normal
data_for_model$prior_overdispersion_mean_association = prior_overdispersion_mean_association$sd_gamma
data_for_model$prior_mean_intercept = prior_mean$intercept_normal
data_for_model$prior_mean_coefficients = prior_mean$coefficients_normal
data_for_model$exclude_priors = exclude_priors
data_for_model$enable_loo = TRUE & enable_loo

Expand All @@ -285,7 +285,12 @@ sccomp_glm_data_frame_counts = function(.data,
verbose = verbose,
seed = mcmc_seed,
max_sampling_iterations = max_sampling_iterations,
pars = c("beta", "alpha", "prec_coeff","prec_sd", "alpha_normalised", "beta_random_intercept", "log_lik")
pars = c("beta", "alpha", "prec_coeff","prec_sd", "alpha_normalised", "beta_random_intercept", "log_lik",

# For getting the priors
"beta_raw_raw"

)
)


Expand Down Expand Up @@ -365,16 +370,93 @@ get_mean_precision = function(fit, data_for_model){
# select(-data, -X)
}

get_mean_precision_association = function(fit){
c(
fit %>%
summary("prec_coeff") %$%
summary %>%
.[,1] ,

fit %>%
summary("prec_sd") %$%
summary %>%
.[,1]
#' Get Mean Precision Association from Fit
#'
#' This function extracts and processes precision-related parameters from a fitted
#' model object (assumed to be from `rstan`). It retrieves the mean and confidence
#' intervals for the intercept and slope priors, as well as the maximum likelihood
#' estimate for the standard deviation of a gamma prior.
#'
#' @param fit A model object, typically from an rstan fit.
#'
#' @return A list containing three elements:
#' - `intercept_normal`: Mean and confidence intervals of the intercept.
#' - `slope_normal`: Mean and confidence intervals of the slope.
#' - `sd_gamma`: Maximum likelihood estimate of the standard deviation
#' of the gamma prior.
#'
#' @importFrom rstan summary
#' @importFrom fitdistrplus fitdist
#' @noRd
get_prior_overdispersion_mean_association <- function(fit){
list(
# Extracting intercept normal prior mean and confidence intervals
intercept_normal =
fit |>
rstan::summary("prec_coeff") %$%
summary |>
_[1, c(1, 3)],

# Extracting slope normal prior mean and confidence intervals
slope_normal =
fit |>
rstan::summary("prec_coeff") %$%
summary |>
_[2, c(1, 3)],

# Fitting a gamma distribution to the extracted standard deviation and obtaining MLE estimates
sd_gamma =
fit |>
rstan::extract("prec_sd") |>
_[[1]] |>
as.numeric() |>
fitdistrplus::fitdist(distr = "gamma", method = "mle") %$%
estimate
)
}

#' Get Mean Precision Association from Fit
#'
#' This function extracts and processes precision-related parameters from a fitted
#' model object (assumed to be from `rstan`). It retrieves the mean and confidence
#' intervals for the intercept and slope priors, as well as the maximum likelihood
#' estimate for the standard deviation of a gamma prior.
#'
#' @param fit A model object, typically from an rstan fit.
#'
#' @return A list containing three elements:
#' - `intercept_normal`: Mean and confidence intervals of the intercept.
#' - `slope_normal`: Mean and confidence intervals of the slope.
#' - `sd_gamma`: Maximum likelihood estimate of the standard deviation
#' of the gamma prior.
#'
#' @importFrom rstan summary
#' @importFrom fitdistrplus fitdist
#' @noRd
get_prior_mean <- function(fit){
list(

# Fitting a gamma distribution to the extracted standard deviation and obtaining MLE estimates
intercept_normal = c(
0,
fit |>
rstan::summary("beta_raw_raw") %$%
summary |>
as_tibble(rownames="parameter") |>
filter(parameter |> str_detect("\\[1")) |>
pull(sd) |>
mean()
),
coefficients_normal = c(
0,
fit |>
rstan::summary("beta_raw_raw") %$%
summary |>
as_tibble(rownames="parameter") |>
filter(parameter |> str_detect("\\[1", negate = TRUE)) |>
pull(sd) |>
mean()
)
)
}
35 changes: 20 additions & 15 deletions R/methods.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,8 +99,8 @@ sccomp_estimate <- function(.data,
bimodal_mean_variability_association = FALSE,
percent_false_positive = 5,
approximate_posterior_inference = "none",
prior_mean = list(intercept = c(0,1), coefficients = c(0,1)),
prior_overdispersion_mean_association = list(intercept = c(5, 2), slope = c(0, 0.6), standard_deviation = c(10, 20)),
prior_mean = list(intercept_normal = c(0,1), coefficients_normal = c(0,1)),
prior_overdispersion_mean_association = list(intercept_normal = c(5, 2), slope_normal = c(0, 0.6), sd_gamma = c(20, 40)),
.sample_cell_group_pairs_to_exclude = NULL,
verbose = TRUE,
enable_loo = FALSE,
Expand All @@ -126,8 +126,8 @@ sccomp_estimate.Seurat = function(.data,
bimodal_mean_variability_association = FALSE,
percent_false_positive = 5,
approximate_posterior_inference = "none",
prior_mean = list(intercept = c(0,1), coefficients = c(0,1)),
prior_overdispersion_mean_association = list(intercept = c(5, 2), slope = c(0, 0.6), standard_deviation = c(10, 20)),
prior_mean = list(intercept_normal = c(0,1), coefficients_normal = c(0,1)),
prior_overdispersion_mean_association = list(intercept_normal = c(5, 2), slope_normal = c(0, 0.6), sd_gamma = c(20, 40)),
.sample_cell_group_pairs_to_exclude = NULL,
verbose = TRUE,
enable_loo = FALSE,
Expand Down Expand Up @@ -186,8 +186,8 @@ sccomp_estimate.SingleCellExperiment = function(.data,
bimodal_mean_variability_association = FALSE,
percent_false_positive = 5,
approximate_posterior_inference = "none",
prior_mean = list(intercept = c(0,1), coefficients = c(0,1)),
prior_overdispersion_mean_association = list(intercept = c(5, 2), slope = c(0, 0.6), standard_deviation = c(10, 20)),
prior_mean = list(intercept_normal = c(0,1), coefficients_normal = c(0,1)),
prior_overdispersion_mean_association = list(intercept_normal = c(5, 2), slope_normal = c(0, 0.6), sd_gamma = c(20, 40)),
.sample_cell_group_pairs_to_exclude = NULL,
verbose = TRUE,
enable_loo = FALSE,
Expand Down Expand Up @@ -249,8 +249,8 @@ sccomp_estimate.DFrame = function(.data,
bimodal_mean_variability_association = FALSE,
percent_false_positive = 5,
approximate_posterior_inference = "none",
prior_mean = list(intercept = c(0,1), coefficients = c(0,1)),
prior_overdispersion_mean_association = list(intercept = c(5, 2), slope = c(0, 0.6), standard_deviation = c(10, 20)),
prior_mean = list(intercept_normal = c(0,1), coefficients_normal = c(0,1)),
prior_overdispersion_mean_association = list(intercept_normal = c(5, 2), slope_normal = c(0, 0.6), sd_gamma = c(20, 40)),
.sample_cell_group_pairs_to_exclude = NULL,
verbose = TRUE,
enable_loo = FALSE,
Expand Down Expand Up @@ -311,8 +311,8 @@ sccomp_estimate.data.frame = function(.data,
bimodal_mean_variability_association = FALSE,
percent_false_positive = 5,
approximate_posterior_inference = "none",
prior_mean = list(intercept = c(0,1), coefficients = c(0,1)),
prior_overdispersion_mean_association = list(intercept = c(5, 2), slope = c(0, 0.6), standard_deviation = c(10, 20)),
prior_mean = list(intercept_normal = c(0,1), coefficients_normal = c(0,1)),
prior_overdispersion_mean_association = list(intercept_normal = c(5, 2), slope_normal = c(0, 0.6), sd_gamma = c(20, 40)),
.sample_cell_group_pairs_to_exclude = NULL,
verbose = TRUE,
enable_loo = FALSE,
Expand Down Expand Up @@ -602,8 +602,6 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate,
)


#fit_model(stan_model("inst/stan/glm_multi_beta_binomial.stan"), chains= 4, output_samples = 500, approximate_posterior_inference = FALSE, verbose = TRUE)

rng2 = rstan::gqs(
stanmodels$glm_multi_beta_binomial_generate_date,
#rstan::stan_model("inst/stan/glm_multi_beta_binomial_generate_date.stan"),
Expand Down Expand Up @@ -693,7 +691,10 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate,
verbose = verbose,
seed = mcmc_seed,
max_sampling_iterations = max_sampling_iterations,
pars = c("beta", "alpha", "prec_coeff","prec_sd", "alpha_normalised", "beta_random_intercept", "log_lik")
pars = c("beta", "alpha", "prec_coeff","prec_sd", "alpha_normalised", "beta_random_intercept", "log_lik",

# For getting the priors
"beta_raw_raw")
)

# Create a dummy tibble
Expand Down Expand Up @@ -856,8 +857,12 @@ sccomp_test.sccomp_tbl = function(.data,

# Add back attributes
add_attr(
.data |> attr("fit") |> get_mean_precision_association(),
"mean_concentration_association"
.data |> attr("fit") |> get_prior_overdispersion_mean_association(),
"prior_overdispersion_mean_association"
) |>
add_attr(
.data |> attr("fit") |> get_prior_mean(),
"prior_mean"
)

if(pass_fit)
Expand Down
6 changes: 4 additions & 2 deletions R/utilities.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -459,7 +459,9 @@ label_deleterious_outliers = function(.my_data){
fit_model = function(
data_for_model, model, censoring_iteration = 1, cores = detectCores(), quantile = 0.95,
warmup_samples = 300, approximate_posterior_inference = TRUE, verbose = FALSE,
seed , pars = c("beta", "alpha", "prec_coeff","prec_sd"), output_samples = NULL, chains=NULL, max_sampling_iterations = 20000
seed , pars = c(
"beta", "alpha", "prec_coeff","prec_sd"
), output_samples = NULL, chains=NULL, max_sampling_iterations = 20000
)
{

Expand Down Expand Up @@ -2123,7 +2125,7 @@ mutate_from_expr_list = function(x, formula_expr, ignore_errors = TRUE){
contrasts_not_in_the_model = contrasts_not_in_the_model[contrasts_not_in_the_model!=""]

if(length(contrasts_not_in_the_model) > 0 & !ignore_errors)
warning(sprintf("sccomp says: These components of your contrasts are not present in the model as parameters: %s", paste(contrasts_not_in_the_model, sep = ", ")))
warning(sprintf("sccomp says: These components of your contrasts are not present in the model as parameters: %s. Factors including special characters, e.g. \"(Intercept)\" require backquotes e.g. \"`(Intercept)`\" ", paste(contrasts_not_in_the_model, sep = ", ")))

# Calculate
if(ignore_errors) my_mutate = mutate_ignore_error
Expand Down
Loading