Skip to content

Commit

Permalink
finished updates to summary function -- #closes #46
Browse files Browse the repository at this point in the history
  • Loading branch information
jr-leary7 committed Sep 20, 2023
1 parent 13d6060 commit c870ab2
Show file tree
Hide file tree
Showing 8 changed files with 240 additions and 186 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: scLANE
Type: Package
Title: Model gene expression dynamics with spline-based NB GLMs, GEEs, & GLMMs
Version: 0.7.2
Version: 0.7.3
Authors@R: c(person(given = "Jack", family = "Leary", email = "j.leary@ufl.edu", role = c("aut", "cre")),
person(given = "Rhonda", family = "Bacher", email = "rbacher@ufl.edu", role = c("ctb", "fnd")))
Description: This package uses truncated power basis spline models to build flexible, interpretable models of single cell gene expression over pseudotime or latent time.
Expand All @@ -10,7 +10,7 @@ Description: This package uses truncated power basis spline models to build flex
License: MIT + file LICENSE
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.2.3
RoxygenNote: 7.2.1
Depends:
glm2,
magrittr,
Expand Down
4 changes: 0 additions & 4 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -39,18 +39,14 @@ importFrom(dplyr,filter)
importFrom(dplyr,group_by)
importFrom(dplyr,if_else)
importFrom(dplyr,inner_join)
importFrom(dplyr,lag)
importFrom(dplyr,lead)
importFrom(dplyr,mutate)
importFrom(dplyr,ntile)
importFrom(dplyr,pull)
importFrom(dplyr,relocate)
importFrom(dplyr,rename)
importFrom(dplyr,rowwise)
importFrom(dplyr,sample_frac)
importFrom(dplyr,select)
importFrom(dplyr,summarise)
importFrom(dplyr,ungroup)
importFrom(dplyr,with_groups)
importFrom(foreach,"%dopar%")
importFrom(foreach,foreach)
Expand Down
20 changes: 16 additions & 4 deletions R/fitGLMM.R
Original file line number Diff line number Diff line change
Expand Up @@ -82,14 +82,20 @@ fitGLMM <- function(X_pred = NULL,
Y = Y,
subject = id.vec,
.before = 1)

if (is.null(Y.offset)) {
glmm_mod <- glmmTMB::glmmTMB(mod_formula,
data = glmm_basis_df,
family = glmmTMB::nbinom2(link = "log"),
se = TRUE,
REML = FALSE)
} else {
glmm_mod <- glmmTMB::glmmTMB(mod_formula,
data = glmm_basis_df,
offset = log(1 / Y.offset),
family = glmmTMB::nbinom2(link = "log"),
se = TRUE,
REML = FALSE)

}
} else {
glmm_basis_df <- data.frame(X1 = tp1(X_pred[, 1], t = round(as.numeric(stats::quantile(X_pred[, 1], 1/3)), 4)),
X2 = tp1(X_pred[, 1], t = round(as.numeric(stats::quantile(X_pred[, 1], 2/3)), 4)),
Expand All @@ -105,14 +111,20 @@ fitGLMM <- function(X_pred = NULL,
paste0("h_PT_", round(as.numeric(stats::quantile(X_pred[, 1], 1/3)), 4)),
paste0("h_", round(as.numeric(stats::quantile(X_pred[, 1], 2/3)), 4), "_PT"),
paste0("h_PT_", round(as.numeric(stats::quantile(X_pred[, 1], 2/3)), 4)))


if (is.null(Y.offset)) {
glmm_mod <- glmmTMB::glmmTMB(Y ~ X1 + X2 + X3 + X4 + (1 + X1 + X2 + X3 + X4 | subject),
data = glmm_basis_df,
family = glmmTMB::nbinom2(link = "log"),
se = TRUE,
REML = FALSE)
} else {
glmm_mod <- glmmTMB::glmmTMB(Y ~ X1 + X2 + X3 + X4 + (1 + X1 + X2 + X3 + X4 | subject),
data = glmm_basis_df,
offset = log(1 / Y.offset),
family = glmmTMB::nbinom2(link = "log"),
se = TRUE,
REML = FALSE)
}
}
# set up results
marge_style_names <- c("B_finalIntercept", marge_style_names)
Expand Down
151 changes: 73 additions & 78 deletions R/summarizeModel.R
Original file line number Diff line number Diff line change
@@ -1,100 +1,95 @@
#' Represent a \code{marge} model as a series of piecewise equations.
#'
#' @name summarizeModel
#' @author Jack Leary & Rhonda Bacher
#' @author Jack Leary
#' @author Rhonda Bacher
#' @import magrittr
#' @importFrom dplyr mutate case_when arrange lead lag rowwise ungroup select
#' @description This function summarizes the model for each gene and allows for quantiative interpretation of fitted gene dynamics..
#' @description This function summarizes the model for each gene and allows for quantiative interpretation of fitted gene dynamics.
#' @param marge.model The fitted model output from \code{\link{marge2}} (this function is internally called by \code{\link{testDynamic}}). Defaults to NULL.
#' @param pseudotime_df The predictor matrix of pseudotime. Defaults to NULL.
#' @param pt The predictor matrix of pseudotime. Defaults to NULL.
#' @return A data.frame of the model coefficients, cutpoint intervals, and formatted equations.
#' @seealso \code{\link{marge2}}
#' @examples
#' \dontrun{
#' summarizeModel(marge.model = marge_mod)
#' summarizeModel(marge.model = marge_mod, pt = pt_df)
#' }

summarizeModel <- function(marge.model = NULL, pt=NULL) {

# check inputs
summarizeModel <- function(marge.model = NULL, pt = NULL) {
# check inputs
if (is.null(marge.model)) { stop("Please provide a non-NULL input argument for marge.model.") }
if (is.null(pt)) { stop("Please provide a non-NULL input argument for pt.") }

if (inherits(marge.model, "try-error")) {
mod_summ <- list(Breakpoint = NA_real_,
mod_summ <- list(Breakpoint = NA_real_,
Slope.Segment = NA_real_,
Trend.Segment = NA_real_)

} else {

# extract model equation & slopes
coef_df <- data.frame(coef_name = names(coef(marge.model$final_mod)),
coef_value = unname(coef(marge.model$final_mod)))

coef_df <- coef_df[-which(coef_df$coef_name == "Intercept"),]

coef_df <- cbind(coef_df, extractBreakpoints(marge.model))

coef_df <- coef_df[order(coef_df$Breakpoint, coef_df$Direction),]
coef_df

MIN <- min(pt[,1])
MAX <- max(pt[,1])
} else {
# extract model equation & slopes
coef_df <- data.frame(coef_name = names(coef(marge.model$final_mod)),
coef_value = unname(coef(marge.model$final_mod)))

coef_ranges <- mapply(function(x,y) {
if (x == "Right") {
c(y, MAX)
} else if (x == "Left") {
c(MIN, y)
}
}, coef_df$Direction, coef_df$Breakpoint)
coef_ranges
# coef_ranges <- coef_ranges[,order(colnames(coef_ranges))]

num_segments <- length(unique(coef_df$Breakpoint)) + 1
mod_seg <- list()
for(i in 1:num_segments) {
if(i == 1) {mod_seg[[i]] <- c(MIN, coef_df$Breakpoint[i])
} else if(i == num_segments) {
mod_seg[[i]] <- c(coef_df$Breakpoint[i-1], MAX)
} else {
mod_seg[[i]] <- c(coef_df$Breakpoint[i-1], coef_df$Breakpoint[i])
}
}
coef_ranges <- t(coef_ranges)
mod_seg_overlaps <- lapply(mod_seg, function(x) {
overlap_ind <- c()
for(i in 1:nrow(coef_ranges)) {
if (x[1] < coef_ranges[i, 2] && x[2] > coef_ranges[i, 1]) {
overlap_ind <- c(overlap_ind, i)
coef_df <- coef_df[-which(coef_df$coef_name == "Intercept"), ]

coef_df <- cbind(coef_df, extractBreakpoints(marge.model))

coef_df <- coef_df[order(coef_df$Breakpoint, coef_df$Direction), ]

MIN <- min(pt[, 1])
MAX <- max(pt[, 1])

coef_ranges <- mapply(function(x, y) {
if (x == "Right") {
c(y, MAX)
} else if (x == "Left") {
c(MIN, y)
}
}, coef_df$Direction, coef_df$Breakpoint)

num_segments <- length(unique(coef_df$Breakpoint)) + 1
mod_seg <- list()
for (i in seq(num_segments)) {
if (i == 1) {
mod_seg[[i]] <- c(MIN, coef_df$Breakpoint[i])
} else if (i == num_segments) {
mod_seg[[i]] <- c(coef_df$Breakpoint[i - 1], MAX)
} else {
mod_seg[[i]] <- c(coef_df$Breakpoint[i - 1], coef_df$Breakpoint[i])
}
}
return(overlap_ind)
})

seg_slopes <- lapply(mod_seg_overlaps, function(x) {

if(length(x) == 1) {
if (coef_df$Direction[x] == "Right") temp_slp <- coef_df$coef_value[x]
if (coef_df$Direction[x] == "Left") temp_slp <- -1*coef_df$coef_value[x]
} else {
to_rev <- which(coef_df$Direction[x] == "Left")
coef_df$coef_value[x][to_rev] <- -1 * coef_df$coef_value[x][to_rev]
temp_slp <- sum(coef_df$coef_value[x])
}
return(temp_slp)
})

seg_slopes <- do.call(c, seg_slopes)

seg_trends <- ifelse(seg_slopes > 0, 1, -1)
seg_trends[seg_slopes==0] <- 0

coef_ranges <- t(coef_ranges)
mod_seg_overlaps <- lapply(mod_seg, function(x) {
overlap_ind <- c()
for (i in seq(nrow(coef_ranges))) {
if (x[1] < coef_ranges[i, 2] && x[2] > coef_ranges[i, 1]) {
overlap_ind <- c(overlap_ind, i)
}
}
return(overlap_ind)
})

seg_slopes <- lapply(mod_seg_overlaps, function(x) {
if (length(x) == 1) {
if (coef_df$Direction[x] == "Right") {
temp_slp <- coef_df$coef_value[x]
} else if (coef_df$Direction[x] == "Left") {
temp_slp <- -1 * coef_df$coef_value[x]
}
} else {
to_rev <- which(coef_df$Direction[x] == "Left")
coef_df$coef_value[x][to_rev] <- -1 * coef_df$coef_value[x][to_rev]
temp_slp <- sum(coef_df$coef_value[x])
}
return(temp_slp)
})
# combine results
seg_slopes <- do.call(c, seg_slopes)
# create discretized trend summary
seg_trends <- ifelse(seg_slopes > 0, 1, -1)
seg_trends[seg_slopes == 0] <- 0
# prepare results
mod_summ <- list(Breakpoint = unique(coef_df$Breakpoint),
Slope.Segment = seg_slopes,
Trend.Segment = seg_trends)

mod_summ <- list(Breakpoint = unique(coef_df$Breakpoint),
Slope.Segment = seg_slopes,
Trend.Segment = seg_trends)
mod_summ

}
return(mod_summ)
}
Loading

0 comments on commit c870ab2

Please sign in to comment.