Added more versions of keystoness. Ranges still off so I tried creating my own MTI function. Made no difference

Author: MGrezlik

.Rhistory

@@ -0,0 +1,35 @@
+ks <- neusRpath::keystoneness("GB")
+ks <- neusRpath::keystoneness("GB")
+neusRpath::plot_keystoneness(ks, top_n = 10)
+export(plot_keystoneness)
+devtools::document()
+devtools::install()
+# To run keystoneness function, EPU required as input
+neusRpath::keystoneness("MAB")
+remove.packages("neusRpath")
+devtools::install()
+# To run keystoneness function, EPU required as input
+neusRpath::keystoneness("MAB")
+# To run plot_keystoneness, need to first run and save keystoneness. Then tell it how many to plot
+ks <- neusRpath::keystoneness("GB")
+neusRpath::plot_keystoneness(ks, top_n = 10)
+range(ks$keystoneness)
+# To run plot_keystoneness, need to first run and save keystoneness. Then tell it how many to plot
+ks <- neusRpath::keystoneness("MAB")
+neusRpath::plot_keystoneness(ks, top_n = 10)
+range(ks$keystoneness)
+# To run plot_keystoneness, need to first run and save keystoneness. Then tell it how many to plot
+ks <- neusRpath::keystoneness("GOM")
+neusRpath::plot_keystoneness(ks, top_n = 10)
+range(ks$keystoneness)
+View(ks)
+devtools::document()
+devtools::install()
+# To run plot_keystoneness, need to first run and save keystoneness. Then tell it how many to plot
+ks <- neusRpath::keystoneness("GB")
+neusRpath::plot_keystoneness(ks, top_n = 10)
+# To run keystoneness function, EPU required as input
+neusRpath::keystoneness("MAB")
+# To run plot_keystoneness, need to first run and save keystoneness. Then tell it how many to plot
+ks <- neusRpath::keystoneness("GB")
+neusRpath::plot_keystoneness(ks, top_n = 10)

NAMESPACE

@@ -1,4 +1,5 @@
 # Generated by roxygen2: do not edit by hand
 
-export(keystoneness)
+export(MTI)
+export(keystoneness_all)
 export(plot_keystoneness)

R/MTI.R

@@ -0,0 +1,67 @@
+#' Calculate Mixed Trophic Impact (MTI)
+#'
+#' @param Rpath Rpath object
+#' @param Rpath.params Rpath parameter object
+#' @param increase Logical; if TRUE, assumes marginal increase
+#'
+#' @return MTI matrix
+#' @export
+MTI <- function(Rpath, Rpath.params, increase = TRUE) {
+  # --- Step 0: prepare ---
+  x <- Rpath.params
+  y <- Rpath
+  
+  all_groups <- x$model$Group
+  n_groups <- length(all_groups)
+  
+  # --- Step 1: Build diet composition (DC) matrix ---
+  DC_full <- matrix(0, nrow = n_groups, ncol = n_groups,
+                    dimnames = list(all_groups, all_groups))
+  
+  preds <- x$diet$Group
+  prey_cols <- colnames(x$diet)[-1]
+  
+  for(i in seq_along(preds)) {
+    predator <- preds[i]
+    if(!predator %in% all_groups) next
+    
+    prey_vals <- as.numeric(dplyr::select(x$diet[i, ], dplyr::all_of(prey_cols)))
+    
+    # keep only prey that exist in all_groups
+    valid_idx <- which(prey_cols %in% all_groups)
+    valid_preys <- prey_cols[valid_idx]
+    valid_vals <- prey_vals[valid_idx]
+    
+    # replace NA with 0
+    valid_vals[is.na(valid_vals)] <- 0
+    
+    DC_full[predator, valid_preys] <- valid_vals
+  }
+  
+  # --- Step 2: Compute fraction of prey consumed by predator (FC) ---
+  bio <- y$BB        # biomass
+  QB  <- y$QB        # consumption/biomass
+  BQB <- bio * QB    # production consumed
+  
+  # Tij = consumption of prey j by predator i
+  Tij <- DC_full * matrix(BQB, nrow = n_groups, ncol = n_groups, byrow = TRUE)
+  
+  Tim <- rowSums(Tij, na.rm = TRUE)  # handle NAs
+  FCij <- matrix(0, nrow = n_groups, ncol = n_groups)
+  
+  for(i in 1:n_groups) {
+    if(!is.na(Tim[i]) && Tim[i] > 0) {
+      FCij[i, ] <- Tij[i, ] / Tim[i]
+    }
+  }
+  
+  FCji <- t(FCij)
+  
+  # --- Step 3: Net impact ---
+  net_impact <- DC_full - FCji
+  
+  # --- Step 4: MTI ---
+  MTI <- MASS::ginv(diag(n_groups) - net_impact) - diag(n_groups)
+  
+  return(MTI)
+}

R/keystoneness.R

@@ -1,70 +1,73 @@
-#' Calculate keystoneness for a Northeast U.S. Rpath model
+#' Calculate keystoneness indices
 #'
-#' Calculates the keystoneness index (KSi) for a given ecosystem model
-#' following Libralato et al. 2006:
-#' KSi = log[epsilon_i * (1 - p_i)],
-#' where epsilon_i is the total impact of group i (from the MTI matrix excluding self-effects)
-#' and p_i is the proportion of total system biomass.
+#' Computes Keystone index #1, #2, #3 and Relative Total Impact (raw and normalized)
+#' for a Northeast U.S. Rpath model, excluding detritus and fleet groups.
 #'
 #' @param epu Character. Ecosystem production unit: "GB", "GOM", or "MAB".
 #' @return A data.frame with ecological groups, biomass, proportion biomass,
-#' total impact, and keystoneness index.
+#' raw and normalized Relative Total Impact, and keystone indices.
 #' @examples
-#' keystoneness("GB")
+#' keystoneness_all("GB")
 #' @export
-keystoneness <- function(epu = c("GB", "GOM", "MAB")) {
+keystoneness_all <- function(epu = c("GB", "GOM", "MAB")) {
 
   epu <- match.arg(epu)
 
   # --- STEP 1: load model and balanced parameters ---
-  data(list = paste0(epu))  # loads GB, GOM, or MAB model into environment
+  data(list = epu)
   model <- get(epu)
 
-  # load balanced parameters
   params_name <- paste0(epu, "_balanced_params")
   data(list = params_name)
   params <- get(params_name)
 
-  # --- STEP 2: Mixed trophic impact matrix ---
-  mti_mat <- Rpath::MTI(Rpath = model, Rpath.params = params)
+  # --- STEP 2: compute MTI matrix using local MTI function ---
+  mti_mat <- MTI(Rpath = model, Rpath.params = params)
 
-  # --- STEP 3: assign group names (ecological groups only) ---
+  # --- STEP 3: filter for living ecological groups only ---
+  eco_rows <- which(params$model$Type < 2)  # living groups only
+  mti_mat <- mti_mat[eco_rows, eco_rows, drop = FALSE]
+  
+  biomass <- model$Biomass[eco_rows]
+  prop_biomass <- biomass / sum(biomass, na.rm = TRUE)
+  
+  # Group names
+  group_names <- model$Group[eco_rows]
   if (!is.null(rownames(model$DC))) {
-    group_names <- rownames(model$DC)
-  } else if (!is.null(model$Group)) {
-    group_names <- model$Group
-  } else {
-    group_names <- paste0("Group", seq_len(nrow(mti_mat)))
+    group_names <- rownames(model$DC)[eco_rows]
   }
 
-  # exclude fishing gear groups
-  n_eco <- nrow(mti_mat) - model$NUM_GEARS
-  mti_mat <- mti_mat[1:n_eco, 1:n_eco]
-  group_names <- group_names[1:n_eco]
+  # --- STEP 4: Relative Total Impact (RTI) ---
+  total_effect <- rowSums(abs(mti_mat), na.rm = TRUE)
+  self_effect  <- diag(abs(mti_mat))
+  rti_raw      <- total_effect - self_effect
+  rti_norm     <- rti_raw / max(rti_raw, na.rm = TRUE)
 
-  # --- STEP 4: compute proportional biomass ---
-  biomass <- model$Biomass[1:n_eco]
-  prop_biomass <- biomass / sum(biomass)
+  # --- STEP 5: Keystone index #1 (Libralato 2006) ---
+  KS1 <- log(rti_raw * (1 - prop_biomass))
 
-  # --- STEP 5: compute total impacts excluding self-effects ---
-  total_effect <- apply(abs(mti_mat), 1, sum)
-  self_effect <- diag(abs(mti_mat))
-  epsilon <- total_effect - self_effect
+  # --- STEP 6: Keystone index #2 (Valls 2015) ---
+  predator_impacts <- rowSums(abs(mti_mat * (mti_mat < 0)), na.rm = TRUE)
+  prey_impacts     <- rowSums(abs(mti_mat * (mti_mat > 0)), na.rm = TRUE)
+  KS2 <- (predator_impacts + prey_impacts) / biomass
 
-  # --- STEP 6: keystoneness (Libralato et al. 2006) ---
-  KS_index <- log(epsilon * (1 - prop_biomass))
+  # --- STEP 7: Keystone index #3 (Valls 2015) ---
+  KS3 <- predator_impacts / biomass
 
-  # --- STEP 7: tidy output ---
+  # --- STEP 8: tidy output ---
   results <- data.frame(
     group = group_names,
     biomass = biomass,
     prop_biomass = prop_biomass,
-    impact = epsilon,
-    keystoneness = KS_index,
+    RTI_raw = rti_raw,
+    RTI_normalized = rti_norm,
+    `Keystone index #1` = KS1,
+    `Keystone index #2` = KS2,
+    `Keystone index #3` = KS3,
     stringsAsFactors = FALSE
   )
 
-  # Attach EPU label as an attribute
+  # Attach EPU label as attribute
   attr(results, "epu") <- switch(epu,
                                  GB  = "Georges Bank",
                                  GOM = "Gulf of Maine",

data-raw/test_runs.R

@@ -5,3 +5,12 @@ neusRpath::keystoneness("MAB")
 ks <- neusRpath::keystoneness("GB")
 neusRpath::plot_keystoneness(ks, top_n = 10)
 
+
+# testing once MTI.R was added
+# From the root of your repo
+source("R/MTI.R")
+source("R/keystoneness.R")
+
+# Test your function
+results <- keystoneness_all("GB")
+head(results)