more functions to glmm (6a10424b) · Commits · github_fork / Presto 20250701

R/glmm.R

+76 −4

Original line number	Diff line number	Diff line
		@@ -52,7 +52,7 @@ glmm_uni <- function(feature, model_base, response, nsim) {


		#' @export
		glmm_multi <- function(formula, design, response, features=NULL, parallel=FALSE, nsim=100, family='poisson') {
		glmm_multi <- function(formula, design, response, features=NULL, parallel=FALSE, nsim=100, family='poisson', return_model=TRUE) {
		if (parallel) {
		future::plan(multicore)
		it_fxn <- furrr::future_map
		@@ -115,6 +115,10 @@ glmm_multi <- function(formula, design, response, features=NULL, parallel=FALSE,
		purrr::reduce(Matrix::rbind2)
		res$response <- response[names(lres), ]

		if (return_model) {
		res$model_base <- model_base
		}

		return(res)
		}

		@@ -162,15 +166,19 @@ compute_I2 <- function(B, S) {


		#' @export
		I2.presto <- function(object, effect, within=NULL) {
		I2.presto <- function(object, effect, effect_levels=NULL, within=NULL, within_levels=NULL) {
		## TODO: save results, check if results were already computed in cache. Option to force recompute.
		if (is.null(within)) {
		if (!effect %in% object$priornames_df$grp) {
		stop(sprintf('Model does not contain random effect intercept for term %s', effect))
		}
		if (is.null(effect_levels)) {
		effect_levels <- subset(object$betanames_df, grpvar == effect)$grp
		}
		idx_use <- object$betanames_df %>%
		tibble::rowid_to_column('idx') %>%
		subset(grpvar == effect) %>%
		# subset(grp %in% effect_levels) %>%
		with(idx)

		I2 <- matrix(compute_I2(t(object$beta[idx_use, ]), t(object$sigma[idx_use, ])), ncol = 1)
		@@ -181,6 +189,9 @@ I2.presto <- function(object, effect, within=NULL) {
		c(paste(within, effect, sep=':'), paste(effect, within, sep=':')),
		object$priornames_df$grp
		)]

		## TODO: add functionalty for subsetting on effect and within levels

		effect_name <- effect_name[!is.na(effect_name)]
		if (length(effect_name) == 0) {
		stop(sprintf('Model does not contain interaction term for %s and %s', effect, within))
		@@ -208,6 +219,7 @@ I2.presto <- function(object, effect, within=NULL) {
		}



		#' @export
		toptable.presto <- function(
		object, n=10, max_pval=.05, max_fdr=.05, min_beta=0, max_beta=Inf, rank_by=c('wald', 'beta', 'pval')[1]
		@@ -251,6 +263,66 @@ toptable.presto <- function(
		}


		#' @export
		genemeans.presto <- function(object, readdepth_varname, xpm=1e6) {
		b0 <- object$beta[which(object$betanames_df$grpvar == '(Intercept)'), ]
		b1 <- object$beta[which(object$betanames_df$grpvar == readdepth_varname), ]
		log_mu <- b0 + log(xpm)*b1
		object$log_mu <- data.frame(log_mu) %>%
		tibble::rownames_to_column('feature')

		return(object)
		}


		#' @export
		effects.presto <- function(object, effect, effect_levels=NULL, within=NULL, within_levels=NULL) {
		## TODO: allow nested effects with the `within` parameter
		## TODO: create user facing accessor function and cache effects computation to not recompute p values

		if (is.null(effect_levels)) {
		effect_levels <- object$betanames_df %>%
		subset(grpvar == effect) %>%
		with(grp)
		}

		object$effects <- dplyr::inner_join(
		cbind(object$betanames_df, object$beta) %>%
		subset(grpvar == effect) %>%
		dplyr::select(-grpvar, -term) %>%
		tidyr::gather(feature, beta, -grp),
		cbind(object$betanames_df, object$sigma) %>%
		subset(grpvar == effect) %>%
		dplyr::select(-grpvar, -term) %>%
		tidyr::gather(feature, sigma, -grp),
		by = c('grp', 'feature')
		) %>%
		dplyr::mutate(
		wald = beta / sigma,
		pval = 2 * (1 - pnorm(abs(beta / sigma)))
		) %>%
		dplyr::mutate(fdr = p.adjust(pval, 'BH')) %>%
		dplyr::left_join(object$log_mu)


		return(object)
		}


		#' @export
		plot_ma.presto <- function(object, max_fdr=.05) {
		object$effect %>%
		ggplot(aes(log_mu, beta, color = fdr < max_fdr)) +
		geom_point(alpha = .8, size = .5) +
		theme_test(base_size = 20) +
		geom_hline(yintercept = 0, linetype = 2, color = 'black') +
		labs(x = 'Expected Mean (logCPM)', y = 'Log Fold Change') +
		scale_color_manual(values = c('grey', 'red')) +
		guides(color = FALSE) +
		NULL

		}


		# #' @export
		# regress_out_one_gene <- function(formula_full, formula_reduced, design, y, common_nUMI) {

Admin message