Dispersion bugfix

NAMESPACE

@@ -1,8 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
-useDynLib(devil, .registration = TRUE)
 export(fit_devil)
 export(test_de)
 importFrom(Rcpp,evalCpp)
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)
+useDynLib(devil);

R/RcppExports.R

@@ -2,46 +2,46 @@
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
 init_beta <- function(y, X) {
-    .Call(`_devil_init_beta`, y, X)
+    .Call('_devil_init_beta', PACKAGE = 'devil', y, X)
 }
 
 beta_fit <- function(y, X, mu_beta, off, max_iter, eps) {
-    .Call(`_devil_beta_fit`, y, X, mu_beta, off, max_iter, eps)
+    .Call('_devil_beta_fit', PACKAGE = 'devil', y, X, mu_beta, off, max_iter, eps)
 }
 
 lte_n_equal_rows <- function(matrix, n, tolerance = 1e-10) {
-    .Call(`_devil_lte_n_equal_rows`, matrix, n, tolerance)
+    .Call('_devil_lte_n_equal_rows', PACKAGE = 'devil', matrix, n, tolerance)
 }
 
 get_row_groups <- function(matrix, n_groups, tolerance = 1e-10) {
-    .Call(`_devil_get_row_groups`, matrix, n_groups, tolerance)
+    .Call('_devil_get_row_groups', PACKAGE = 'devil', matrix, n_groups, tolerance)
 }
 
 make_table_if_small <- function(x, stop_if_larger) {
-    .Call(`_devil_make_table_if_small`, x, stop_if_larger)
+    .Call('_devil_make_table_if_small', PACKAGE = 'devil', x, stop_if_larger)
 }
 
 conventional_loglikelihood_fast <- function(y, mu, log_theta, model_matrix, do_cr_adj, unique_counts = as.numeric( c()), count_frequencies = as.numeric( c())) {
-    .Call(`_devil_conventional_loglikelihood_fast`, y, mu, log_theta, model_matrix, do_cr_adj, unique_counts, count_frequencies)
+    .Call('_devil_conventional_loglikelihood_fast', PACKAGE = 'devil', y, mu, log_theta, model_matrix, do_cr_adj, unique_counts, count_frequencies)
 }
 
 conventional_score_function_fast <- function(y, mu, log_theta, model_matrix, do_cr_adj, unique_counts = as.numeric( c()), count_frequencies = as.numeric( c())) {
-    .Call(`_devil_conventional_score_function_fast`, y, mu, log_theta, model_matrix, do_cr_adj, unique_counts, count_frequencies)
+    .Call('_devil_conventional_score_function_fast', PACKAGE = 'devil', y, mu, log_theta, model_matrix, do_cr_adj, unique_counts, count_frequencies)
 }
 
 conventional_deriv_score_function_fast <- function(y, mu, log_theta, model_matrix, do_cr_adj, unique_counts = as.numeric( c()), count_frequencies = as.numeric( c())) {
-    .Call(`_devil_conventional_deriv_score_function_fast`, y, mu, log_theta, model_matrix, do_cr_adj, unique_counts, count_frequencies)
+    .Call('_devil_conventional_deriv_score_function_fast', PACKAGE = 'devil', y, mu, log_theta, model_matrix, do_cr_adj, unique_counts, count_frequencies)
 }
 
 compute_hessian <- function(beta, overdispersion, y, design_matrix, size_factors) {
-    .Call(`_devil_compute_hessian`, beta, overdispersion, y, design_matrix, size_factors)
+    .Call('_devil_compute_hessian', PACKAGE = 'devil', beta, overdispersion, y, design_matrix, size_factors)
 }
 
 compute_scores <- function(design_matrix, y, beta, overdispersion, size_factors) {
-    .Call(`_devil_compute_scores`, design_matrix, y, beta, overdispersion, size_factors)
+    .Call('_devil_compute_scores', PACKAGE = 'devil', design_matrix, y, beta, overdispersion, size_factors)
 }
 
 compute_clustered_meat <- function(design_matrix, y, beta, overdispersion, size_factors, clusters) {
-    .Call(`_devil_compute_clustered_meat`, design_matrix, y, beta, overdispersion, size_factors, clusters)
+    .Call('_devil_compute_clustered_meat', PACKAGE = 'devil', design_matrix, y, beta, overdispersion, size_factors, clusters)
 }
 

R/main.R

@@ -21,6 +21,7 @@
 #' @return List containing fitted model parameters including beta coefficients,
 #' overdispersion parameter (if estimated), and beta sigma.
 #' @export
+#' @rawNamespace useDynLib(devil);
 fit_devil <- function(input_matrix, design_matrix, overdispersion = TRUE, offset=0, size_factors=TRUE, verbose=FALSE, max_iter=500, eps=1e-4) {
 
   if (size_factors) {
@@ -77,7 +78,7 @@ fit_devil <- function(input_matrix, design_matrix, overdispersion = TRUE, offset
     #   fit_dispersion(beta[i,], design_matrix, input_matrix[i,])
     # }) %>% unlist()
   } else {
-    theta = 0
+    theta = rep(0, ngenes)
   }
 
   return(list(beta=beta, overdispersion=theta, sigma_beta=sigma, iterations=iterations, size_factors=sf, offset_matrix=offset_matrix, design_matrix=design_matrix, input_matrix=input_matrix))

R/overdispersion.R

@@ -10,6 +10,9 @@ fit_dispersion <- function(beta, model_matrix, y, offset_matrix, do_cox_reid_adj
   mean_vector[mean_vector == 0] <- 1e-6
   mu <- mean(y)
   start_value <- (stats::var(y) - mu) / mu^2
+  if(is.na(start_value) || start_value <= 0){
+    start_value <- 0.5
+  }
 
   far_left_value <- conventional_score_function_fast(y, mu = mean_vector, log_theta = log(1e-8),
                                                      model_matrix = model_matrix, do_cr_adj = do_cox_reid_adjustment, tab[[1]], tab[[2]])

R/test.R

@@ -11,6 +11,7 @@
 #' @return A tibble containing the results of the differential expression testing.
 #' @details This function computes log fold changes and p-values for each gene in parallel and filters the results based on the maximum absolute log fold change specified.
 #' @export
+#' @rawNamespace useDynLib(devil);
 test_de <- function(devil.fit, contrast, pval_adjust_method = "BH", max_lfc = 10, clusters = NULL) {
   # Extract necessary information
   ngenes <- nrow(devil.fit$input_matrix)
@@ -19,7 +20,8 @@ test_de <- function(devil.fit, contrast, pval_adjust_method = "BH", max_lfc = 10
   # Calculate log fold changes
   lfcs <- (devil.fit$beta %*% contrast) %>%
     unlist() %>%
-    unname()
+    unname() %>%
+    c()
 
   # Calculate p-values in parallel
   if (!is.null(clusters)) {