Correctly group linear modeling

qhtcp-workflow/apps/r/calculate_interaction_zscores.R

@@ -825,22 +825,31 @@ generate_rank_plot_configs <- function(df_filtered, variables, is_lm = FALSE) {
         rectangles <- NULL
       }
       
-      # Check if there is sufficient variation
+      # Fit linear model grouped by OrfRep, Gene, num
-      if (length(unique(df_filtered[[x_var]])) < 2 || length(unique(df_filtered[[y_var]])) < 2) {
+      lm_results <- df_filtered %>%
-        print(df_filtered %>% select(all_of(c("scan", "Plate", "Row", "Col", "num", "conc_num", x_var, y_var))))
+        group_by(OrfRep, Gene, num) %>%
-        message("Not enough variation in ", x_var, " or ", y_var, " to fit linear model. Skipping.")
+        do({
-        next
+          model <- try(lm(as.formula(paste(y_var, "~", x_var)), data = .), silent = TRUE)
+          if (inherits(model, "try-error")) {
+            # Return NA if model fails
+            data.frame(intercept = NA, slope = NA, r_squared = NA)
+          } else {
+            summary_model <- summary(model)
+            data.frame(
+              intercept = coef(model)[1],
+              slope = coef(model)[2],
+              r_squared = summary_model$r.squared
+            )
           }
+        }) %>%
+        ungroup()
       
-      # Fit linear model
+      aggregated_lm <- lm_results %>%
-      lm_fit <- lm(df_filtered[[y_var]] ~ df_filtered[[x_var]], data = df_filtered)
+        summarize(
-      
+          intercept = mean(intercept, na.rm = TRUE),
-      # # Check for perfect fit
+          slope = mean(slope, na.rm = TRUE),
-      # if (summary(lm_fit)$sigma == 0) {
+          r_squared = mean(r_squared, na.rm = TRUE)
-      #   next  # Skip this iteration if the fit is perfect
+        )
-      # }
-      
-      r_squared <- summary(lm_fit)$r.squared
       
       configs[[length(configs) + 1]] <- list(
         df = df_filtered,
@@ -852,14 +861,14 @@ generate_rank_plot_configs <- function(df_filtered, variables, is_lm = FALSE) {
           list(
             x = 0,
             y = 0,
-            label = paste("R-squared =", round(r_squared, 2))
+            label = paste("R-squared =", round(aggregated_lm$r_squared, 2))
           )
         ),
         sd_band_values = NULL, # Not applicable
         shape = 3,
         size = 0.1,
         add_smooth = TRUE,
-        lm_line = list(intercept = coef(lm_fit)[1], slope = coef(lm_fit)[2]),
+        lm_line = list(intercept = aggregated_lm$intercept, slope = aggregated_lm$slope),
         legend_position = "right",
         color_var = "Overlap",
         x_label = x_var,
@@ -1273,14 +1282,14 @@ main <- function() {
         ungroup()
 
       message("Calculating reference strain interaction scores")
-      reference_results <- calculate_interaction_scores(reference_strain, max_conc, interaction_vars, group_vars = c("OrfRep", "Gene", "num"))
+      reference_results <- calculate_interaction_scores(reference_strain, max_conc, group_vars = c("OrfRep", "Gene", "num"))
       zscores_calculations_reference <- reference_results$calculations
       zscores_interactions_reference <- reference_results$interactions
       zscores_calculations_reference_joined <- reference_results$calculations_joined
       zscores_interactions_reference_joined <- reference_results$interactions_joined
 
       message("Calculating deletion strain(s) interactions scores")
-      deletion_results <- calculate_interaction_scores(deletion_strains, max_conc, interaction_vars, group_vars = c("OrfRep", "Gene", "num"))
+      deletion_results <- calculate_interaction_scores(deletion_strains, max_conc, group_vars = c("OrfRep", "Gene", "num"))
       zscores_calculations <- deletion_results$calculations
       zscores_interactions <- deletion_results$interactions
       zscores_calculations_joined <- deletion_results$calculations_joined