.Rhistory

} else {
group_names <- paste0("Group", seq_len(nrow(mti_mat)))
}
# 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: compute proportional biomass ---
biomass <- model$Biomass[1:n_eco]
prop_biomass <- biomass / sum(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: keystoneness (Libralato et al. 2006) ---
KS_index <- log(epsilon * (1 - prop_biomass))
KS_index
range(KS_index)
max(KS_index)
ks <- neusRpath::keystoneness("GB")
ks <- neusRpath::keystoneness("GB")
neusRpath::plot_keystoneness(ks, top_n = 10)
# 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)
load("C:/Users/maxwell.grezlik/Downloads/SurvdatBio_EPU.RData")
View(survdat.bio.epu)
Survdat.bio <- readRDS("C:/Users/maxwell.grezlik/Downloads/Survdat.bio.rds")
Survdat_bio_EPU <- readRDS("C:/Users/maxwell.grezlik/Downloads/Survdat_bio_EPU.rds")
View(Survdat_bio_EPU)
colnames(Survdat_bio_EPU)
load(here::here('data/GB.rda'))
load(here::here('data/GB_balanced_params.rda'))
MTI_out <- MTI_new(GB, GB_balanced_params)
# MTG 01/27/2026
# keystoneness values calculated as described in the literature do not align with values
# calculated in Ecopath.
# Kerim found old code from Villy which corrects the MTI function
# MTI_new function is copied from his rpath_notes repo
# That repo is public for now but will need to be private
library(Rpath)
MTI_new <- function(Rpath, Rpath.params){
# Useful indices
N_BIO <- Rpath$NUM_LIVING+Rpath$NUM_DEAD
# Make empty matrices
AA <- matrix(0, Rpath$NUM_GROUPS, Rpath$NUM_GROUPS)
colnames(AA) <- Rpath$Group
rownames(AA) <- Rpath$Group
DCC <- AA
# Diet matrix removing import row, group column, adding detrital columns
DCC[1:N_BIO,1:N_BIO] <- cbind(as.matrix(Rpath.params$diet[,-"Group"])[1:(dim(Rpath.params$diet)[1]-1),],
matrix(0, nrow=dim(Rpath.params$diet)[1]-1, ncol=Rpath$NUM_DEAD))
# Landings only not discards in DCC (discards never included at all in this??)
DCC[,(N_BIO+1):Rpath$NUM_GROUPS] <- Rpath$Landings
# Normalizing and cleaning DCC (diet composition and landings)
DCC[is.na(DCC)] <- 0
DCC <- t(t(DCC)/colSums(DCC))
DCC[is.nan(DCC) | is.infinite(DCC)] <- 0
# Living total flows (and flows FROM detritus via diets, but not TO detritus)
for (i in 1:N_BIO){
for (j in 1:Rpath$NUM_LIVING){
AA[i,j] <- Rpath$Biomass[j] * Rpath$QB[j] * DCC[i,j]
}
}
# These lines not needed because flow into detritus in AA matrix is 0'd out below?
# Keeping in comments for reference.
#  detflow <- (Rpath$PB*Rpath$Biomass*(1-Rpath$EE)) + Rpath$QB*Rpath$Biomass*Rpath$Unassim
#  detmat <- detflow * Rpath$DetFate
#  AA[,(Rpath$NUM_LIVING+1):N_BIO] <- detmat
# Should Discards be here?  Seems to work but check later?
AA[,(N_BIO+1):Rpath$NUM_GROUPS] <- Rpath$Landings + Rpath$Discards
fishing_det <- colSums(Rpath$Discards) * Rpath$DetFate[(N_BIO+1):Rpath$NUM_GROUPS,]
AA[(N_BIO+1):Rpath$NUM_GROUPS, (Rpath$NUM_LIVING+1):N_BIO] <- fishing_det
# "Detritus is not considered a predator" - remove flows to detritus before normalizing
AA[,(Rpath$NUM_LIVING+1):N_BIO] <- 0
# Normalize FC and and then flip so that predators are in columns
FC <- AA/rowSums(AA)
FC[is.nan(FC) | is.infinite(FC)] <- 0
FC_T <- t(FC)
# Direct impacts matrix
MTI <- (DCC - FC_T)
# From Ulanowicz paper, find inverse of ([Identity] - MTI)
MTI_diag <- diag(Rpath$NUM_GROUPS) - MTI
# Invert the matrix
MTI_inv <- MASS::ginv((MTI_diag))
# Now subtract 1 off the diagonal matrix
diag(MTI_inv) <- diag(MTI_inv) - 1
# The final(?) piece, hidden in Villy's '03 code - 0 out the detritial Identity
for (i in (Rpath$NUM_LIVING+1):N_BIO){
MTI_inv[i,i] <- 0
}
colnames(MTI_inv) <- Rpath$Group
rownames(MTI_inv) <- Rpath$Group
return(MTI_inv)
}
MTI_out <- MTI_new(GB, GB_balanced_params)
#'
#' 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,
#' raw and normalized Relative Total Impact, and keystone indices.
#' @examples
#' keystoneness_all("GB")
#' @export
keystoneness_all <- function(epu = c("GB", "GOM", "MAB")) {
epu <- match.arg(epu)
# --- STEP 1: load model and balanced parameters ---
data(list = epu)
model <- get(epu)
params_name <- paste0(epu, "_balanced_params")
data(list = params_name)
params <- get(params_name)
# --- STEP 2: compute MTI matrix using local MTI function ---
mti_mat <- MTI(Rpath = model, Rpath.params = params)
# --- 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)[eco_rows]
}
# --- 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 5: Keystone index #1 (Libralato 2006) ---
KS1 <- log(rti_raw * (1 - prop_biomass))
# --- 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 7: Keystone index #3 (Valls 2015) ---
KS3 <- predator_impacts / biomass
# --- STEP 8: tidy output ---
results <- data.frame(
group = group_names,
biomass = biomass,
prop_biomass = prop_biomass,
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 attribute
attr(results, "epu") <- switch(epu,
GB  = "Georges Bank",
GOM = "Gulf of Maine",
MAB = "Mid Atlantic Bight"
)
return(results)
}
keystoneness_all(epu = GB)
keystoneness_all(epu = 'GB')
Rpath::MTI(GB)
Rpath::MTI(GB, GB_balanced_params)
MIT_out
MTI_out
MTI_Rpath <- Rpath::MTI(GB, GB_balanced_params)
final_diff <- MTI_Rpath - MTI_out
max(final_diff)
# From Valls et al
epsilon_mat <- MTI_out[1:bal$NUM_LIVING, 1:bal$NUM_LIVING]
GB$NUM_LIVING
# From Valls et al
epsilon_mat <- MTI_out[1:GB$NUM_LIVING, 1:GB$NUM_LIVING]
diag(epsilon_mat) <- 0
epsilon_i <- sqrt(rowSums(epsilon_mat * epsilon_mat))
p_i <- bal$Biomass[1:GB$NUM_LIVING]/sum(GB$Biomass[1:GB$NUM_LIVING])
p_i <- GB$Biomass[1:GB$NUM_LIVING]/sum(GB$Biomass[1:GB$NUM_LIVING])
KS_1 <- log10(epsilon_i * (1-p_i))
KS_2 <- log10(epsilon_i * (1/p_i))
# Biomass Rank in descending order
drank <- GB$NUM_LIVING - rank(GB$Biomass[1:GB$NUM_LIVING]) + 1
KS_3 <- log10(epsilon_i * drank)
Rel.Tot <- (rowSums((MTI_out[1:GB$NUM_LIVING, 1:GB$NUM_LIVING])))
KS_out <- cbind(KS_1,KS_2,KS_3)
KS_out
epu = 'GB'
# --- STEP 1: load model and balanced parameters ---
data(list = epu)
model <- get(epu)
params_name <- paste0(epu, "_balanced_params")
data(list = params_name)
params <- get(params_name)
# --- STEP 2: compute MTI matrix using local MTI function ---
source(here::here('R/mti_new.R'))
mti_mat <- MTI_new(Rpath = model, Rpath.params = params)
#'
#' 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,
#' raw and normalized Relative Total Impact, and keystone indices.
#' @examples
#' keystoneness_all("GB")
#' @export
keystoneness_all <- function(epu = c("GB", "GOM", "MAB")) {
epu <- match.arg(epu)
# --- STEP 1: load model and balanced parameters ---
data(list = epu)
model <- get(epu)
params_name <- paste0(epu, "_balanced_params")
data(list = params_name)
params <- get(params_name)
# --- STEP 2: compute MTI matrix using local MTI function ---
source(here::here('R/mti_new.R'))
mti_mat <- MTI_new(Rpath = model, Rpath.params = params)
# --- 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)[eco_rows]
}
# --- 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 5: Keystone index #1 (Libralato 2006) ---
KS1 <- log(rti_raw * (1 - prop_biomass))
# --- 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 7: Keystone index #3 (Valls 2015) ---
KS3 <- predator_impacts / biomass
# --- STEP 8: tidy output ---
results <- data.frame(
group = group_names,
biomass = biomass,
prop_biomass = prop_biomass,
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 attribute
attr(results, "epu") <- switch(epu,
GB  = "Georges Bank",
GOM = "Gulf of Maine",
MAB = "Mid Atlantic Bight"
)
return(results)
}
keystoneness_all('GB')
keystoneness_all('GOM')
keystoneness_all('MAB')
devtools::load_all()
devtools::load_all()
devtools::load_all()
keystoneness_all('GB')
devtools::document()
#'
#' 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,
#' raw and normalized Relative Total Impact, and keystone indices.
#' @examples
#' keystoneness_all("GB")
#' @export
keystoneness_all <- function(epu = c("GB", "GOM", "MAB")) {
epu <- match.arg(epu)
# --- STEP 1: load model and balanced parameters ---
data(list = epu)
model <- get(epu)
params_name <- paste0(epu, "_balanced_params")
data(list = params_name)
params <- get(params_name)
# --- STEP 2: compute MTI matrix using local MTI function ---
mti_mat <- MTI_new(Rpath = model, Rpath.params = params)
# --- 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)[eco_rows]
}
# --- 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 5: Keystone index #1 (Libralato 2006) ---
KS1 <- log(rti_raw * (1 - prop_biomass))
# --- 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 7: Keystone index #3 (Valls 2015) ---
KS3 <- predator_impacts / biomass
# --- STEP 8: tidy output ---
results <- data.frame(
group = group_names,
biomass = biomass,
prop_biomass = prop_biomass,
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 attribute
attr(results, "epu") <- switch(epu,
GB  = "Georges Bank",
GOM = "Gulf of Maine",
MAB = "Mid Atlantic Bight"
)
return(results)
}
keystoneness_all('GB')
#' @param epu Character. Ecosystem production unit:
#'   one of "GB", "GOM", or "MAB".
#'
#' @return A data.frame with ecological group names, biomass,
#' proportion of biomass, raw and normalized RTI, and keystone indices.
#'
#' @examples
#' keystoneness_all("GB")
#'
#' @export
keystoneness_all <- function(epu = c("GB", "GOM", "MAB")) {
epu <- match.arg(epu)
## ------------------------------------------------------------------
## STEP 1: Load balanced model and parameters
## ------------------------------------------------------------------
# Load model object (GB, GOM, or MAB)
data(list = epu, envir = environment())
model <- get(epu, envir = environment())
# Load corresponding balanced parameters
params_name <- paste0(epu, "_balanced_params")
data(list = params_name, envir = environment())
params <- get(params_name, envir = environment())
## ------------------------------------------------------------------
## STEP 2: Compute MTI matrix (Ecopath-consistent)
## ------------------------------------------------------------------
# MTI_new() is assumed to live in R/mti_new.R and be loaded with the package
mti_mat <- MTI_new(
Rpath        = model,
Rpath.params = params
)
## ------------------------------------------------------------------
## STEP 3: Restrict to living ecological groups
## ------------------------------------------------------------------
eco_rows <- which(params$model$Type < 2)
# Defensive check: living groups should match NUM_LIVING
if (length(eco_rows) != model$NUM_LIVING) {
warning(
"Number of living groups (Type < 2) does not match model$NUM_LIVING"
)
}
mti_mat <- mti_mat[eco_rows, eco_rows, drop = FALSE]
biomass <- model$Biomass[eco_rows]
group_names <- model$Group[eco_rows]
prop_biomass <- biomass / sum(biomass, na.rm = TRUE)
## ------------------------------------------------------------------
## STEP 4: Relative Total Impact (RTI)
## ------------------------------------------------------------------
total_effect <- rowSums(abs(mti_mat), na.rm = TRUE)
self_effect  <- abs(diag(mti_mat))
rti_raw <- total_effect - self_effect
# Normalize by maximum (explicit choice)
rti_norm <- rti_raw / max(rti_raw, na.rm = TRUE)
## ------------------------------------------------------------------
## STEP 5: Keystone index #1 (Libralato et al. 2006)
## ------------------------------------------------------------------
# Guard against zeros and log(0)
eps <- .Machine$double.eps
KS1 <- log10(pmax(rti_raw * (1 - prop_biomass), eps))
## ------------------------------------------------------------------
## STEP 6: Keystone indices #2 and #3 (Valls et al. 2015-style partitioning)
## ------------------------------------------------------------------
predator_impacts <- rowSums(abs(mti_mat * (mti_mat < 0)), na.rm = TRUE)
prey_impacts     <- rowSums(abs(mti_mat * (mti_mat > 0)), na.rm = TRUE)
# Avoid division by zero
biomass_safe <- pmax(biomass, eps)
KS2 <- (predator_impacts + prey_impacts) / biomass_safe
KS3 <- predator_impacts / biomass_safe
## ------------------------------------------------------------------
## STEP 7: Tidy output
## ------------------------------------------------------------------
results <- data.frame(
group             = group_names,
biomass           = biomass,
prop_biomass      = prop_biomass,
RTI_raw           = rti_raw,
RTI_normalized    = rti_norm,
Keystone_index_1  = KS1,
Keystone_index_2  = KS2,
Keystone_index_3  = KS3,
stringsAsFactors  = FALSE
)
attr(results, "epu") <- switch(
epu,
GB  = "Georges Bank",
GOM = "Gulf of Maine",
MAB = "Mid Atlantic Bight"
)
return(results)
}
keystoneness_all('GB')
devtools::load_all()
keystoneness_all('GB')
summary(keystoneness_all("GB"))
devtools::load_all()
keystoneness_all('GB')
summary(keystoneness_all("GB"))
devtools::load_all()
keystoneness_all('GB')
keystoneness_all('GB')
devtools::load_all()
keystoneness_all('GB')
library(neusRpath)
out_dir <- file.path("inst", "extdata", "keystoneness")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
for (epu in c("GB", "GOM", "MAB")) {
ks <- keystoneness_all(epu)
write.csv(
ks,
file = file.path(out_dir, paste0(epu, "_keystoneness.csv")),
row.names = FALSE
)
}
source("data-raw/write_keystoneness_csvs.R")
plot_keystoneness('GB')
# 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")
# To run plot_keystoneness, need to first run and save keystoneness. Then tell it how many to plot
ks <- neusRpath::keystoneness_all("GB")
neusRpath::plot_keystoneness(ks, top_n = 10)
# Test your function
results <- keystoneness_all("GB")
head(results)
library(readr)
GB_keystoneness <- read_csv("inst/extdata/keystoneness/GB_keystoneness.csv")
View(GB_keystoneness)
summary(GB_keystoneness)
devtools::load_all()
devtools::document()
grep -R "plot_keystoneness" .
devtools::document()
devtools::load_all()
library(neusRpath)
out_dir <- file.path("inst", "extdata", "keystoneness")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
for (epu in c("GB", "GOM", "MAB")) {
ks <- keystoneness_all(epu)
write.csv(
ks,
file = file.path(out_dir, paste0(epu, "_keystoneness.csv")),
row.names = FALSE
)
}
library(neusRpath)
out_dir <- file.path("inst", "extdata", "keystoneness")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
for (epu in c("GB", "GOM", "MAB")) {
ks <- keystoneness_all(epu)
write.csv(
ks,
file = file.path(out_dir, paste0(epu, "_keystoneness.csv")),
row.names = FALSE
)
}
devtools::load_all()
library(neusRpath)
out_dir <- file.path("inst", "extdata", "keystoneness")
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
for (epu in c("GB", "GOM", "MAB")) {
ks <- keystoneness_all(epu)
write.csv(
ks,
file = file.path(out_dir, paste0(epu, "_keystoneness.csv")),
row.names = FALSE
)
}
library(readr)
GB_keystoneness <- read_csv("inst/extdata/keystoneness/GB_keystoneness.csv")
View(GB_keystoneness)
summary(GB_keystoneness)
library(readr)
GOM_keystoneness <- read_csv("inst/extdata/keystoneness/GOM_keystoneness.csv")
View(GOM_keystoneness)
library(readr)
MAB_keystoneness <- read_csv("inst/extdata/keystoneness/MAB_keystoneness.csv")
View(MAB_keystoneness)
summary(GOM_keystoneness)
summary(GB_keystoneness)
summary(MAB_keystoneness)
.Machine$double.eps