In NOAA-EDAB/RpathQNM: What the Package Does (One Line, Title Case)
Setup
require(devtools)
require(Rtools)
require(data.table)
require (Rpath)
require(QPress)
require(tidyverse)
# install_github("NOAA-EDAB/Rpath")
# install_github("SWotherspoon/QPress")
Rpath for WSS28
#Western Scotian Shelf Model (No stanza groups)
#SML
#User parameters
#-------------------------------------------------------------------------------
#Required packages
library(data.table); library(Rpath)
setwd("~/JCT work/WGNARS/RpathtoQNM/RpathQNM/WSS28model")
#-------------------------------------------------------------------------------
#User created functions
#-------------------------------------------------------------------------------
#Western Scotian Shelf Groups - collapsed list (no multistanza groups)
groups <- c('Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks',
'Large pelagics',
'D piscivores', 'L benthivores', 'Skates', 'Small pelagics', 'Other pelagic mesopelagic',
'Squids', 'Megabenthos', 'Small crab other arthropoda', 'Shrimps',
'Bivalves', 'Other molluscs', 'Worms', 'Meiofauna','Gelatinous zoop',
'Macrozoop', 'Mesozoop', 'Microzoop', 'Microflora', 'Phytoplankton',
'Discards', 'Detritus', 'Fishery')
#Identify type (0 = consumer, 1 = producer, 2 = detritus, 3 = fleet)
types <- c(rep(0, 24), 1, 2, 2, 3)
#Create parameter list object
WSS28.params <- create.rpath.params(groups, types)
#Fill in model parameters
#Biomass, Production, consumption
biomass <- c(0.407075, 0.049754, 0.04383, 0.00617, 0.0260249,0.023573, 4.082286, 4.240253,
0.137757, 6.294171, 0.878285, 0.173893,
0.458598,
1.818481,
1.302563,
64.01978,
2,
7.345485,
42.75382,
0.5204,
41.31,
23.79555,
5.800508,
3.678861,
33.664,
0.063293,
1,NA)
pb <- c(0.071,0.18,0.147007,0.25,0.18,0.4,0.5341109,0.3798308,0.2320817,0.6959251,0.7432774,4,
0.8221899,2.759435,3,0.69,0.75,1.25,1.32583,15.51,3.04,29.2,
82.8, 104.938, 70.639, rep(NA, 3))
qb <- c(4.94,14.5,7.338632,87.6,4.78,4.24,2.868397,3.254083,2.447539,3.475426,3.532147,11.33333,
5.481265,18.39623 ,rep(NA, 4), 8.83887,62.05, 19.5, 73, rep(NA, 6))
WSS28.params$model[, Biomass := biomass]
WSS28.params$model[, PB := pb]
WSS28.params$model[, QB := qb]
#Production to Consumption
WSS28.params$model[Group %in% c('megabenthos', 'Shrimps',
'Bivalves', 'Other molluscs', 'Worms' ),
ProdCons := 0.15]
#WSS28.params$model[Group == 'Mesozoop', ProdCons := 0.4]
WSS28.params$model[Group %in% c('Microzoop', 'Microflora'), ProdCons := 0.5]
#Biomass accumulation and unassimilated production
ba <- c(0.0032566,0,0.003702101,0,-0.001951868,0,-0.05736,0.044246, rep(0, 4), 0.0301295, 0,0,
0.03234285,rep(0,11), NA)
# ba<-c(rep(0,28), NA)
WSS28.params$model[, BioAcc := ba]
WSS28.params$model[, Unassim := c(rep(0.2, 20), rep(0.4, 3), 0.2, rep(0, 3), NA)]
#Detrital fate
WSS28.params$model[, Detritus := c(rep(1, 25), rep(0, 3))]
WSS28.params$model[, Discards := c(rep(0, 27), 1)]
#Landings/Discards
land <- c(rep(0, 4),0.00209077,
0.005067847,
0.3251713,
0.182276,
0.000799,
0.999061,
0.019708,
0.001637,
0.220424,0,
0.000162,
0.117151,
0.00236,0,
0.03026,rep(0,8),NA)
disc <- c(0,0.002995,0,2.4E-05,
0.000575,
0.000469,
0.013477,
0.002266,
0.013094,
0.000487,
0.00056,
1.9E-05,
0.008776,
0.000633,
8E-06 ,0,4.1E-05,0,0.019869, rep(0, 8), NA)
WSS28.params$model[, Fishery := land]
WSS28.params$model[, Fishery.disc := disc]
#Diet composition
DC.groups <- groups[which(!groups %in% c('Discards', 'Detritus', 'Fishery'))]
# DC.data <- as.data.table(matrix(c(0, 0, 0, 0 ,0, 0 ,0, 0, 0, 0, 0,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0, 0,0,0,0,0,0,
# 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0,
# 0, 0 ,0, 0 ,0.006511993, 0 ,0, 0 ,0, 0 ,0, 0, 0, 0 ,0, 0 ,0, 0, 0, 0 ,0, 0, 0, 0, 0,
# 0, 0, 0, 0.002663992, 0, 0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0 ,0, 0 ,0, 0 ,0, 0, 0,
# 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0,
# 0, 0, 0, 0 ,0, 0.002219998 ,0, 0, 0, 0, 0, 0 ,0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0, 0, 0, 0,
# 0.004802, 0.2437588, 0.1276689, 0.0007209978, 0.1848318, 0.06917194, 0.05424953, 0.005541045 ,0.003080484, 0 ,0, 0.003657, 0, 0, 0, 0 ,0 ,0, 0, 0,0, 0, 0, 0 ,0,
# 0.012545, 0.09793591, 0.2855617, 0.005076986, 0.3280216, 0.1511799, 0.04674805, 0.004132936, 0.05460171, 0, 0, 0.007678, 0.001092015, 0, 0, 0 ,0, 0 ,0 ,0, 0, 0 ,0, 0, 0,
# 0, 0 ,0.006431994, 0, 0.0008079991, 0.002838997, 0.0004965155, 0 ,0, 0 ,0, 0.000187, 0, 0 ,0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0,
# 0.173918, 0.3776346 ,0.5207465, 0.3529089, 0.4223085, 0.6149244, 0.1726697, 0.00143368, 0.02304704, 0, 0, 0.015327, 0, 0, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0, 0, 0,
# 0.0194, 0.007396993, 0.05680194, 0.07901775, 0.04387695, 0.1276559, 0.02289399, 0.01478033, 0.03303853, 0, 0, 0.004734, 0 ,0, 0, 0 ,0 ,0 ,0, 0, 0 ,0, 0 ,0, 0,
# 0.011769, 0.2679567, 0.002056998, 0.03595289, 0.01354798, 0.01555398, 0.02838177, 0.001817405 ,0.006008342, 0, 0, 0.037929, 0.002280137, 0, 0, 0, 0, 0 ,0, 0, 0, 0 ,0, 0, 0,
# 0, 0 ,2.099998E-05, 0, 9.29999E-05, 7.299994E-05, 0.003771638, 0.001804974, 0.00659312, 0, 0, 0, 0.007052792, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
# 0, 0, 0, 0, 0, 0, 0.03455577, 0.07374978, 0.1577646, 0.02415398, 0.060237, 0.014737, 0.09158555, 0.009690308, 0.004533, 0, 0 ,0.03769, 0.02, 0.0001840747, 0, 0, 0, 0, 0,
# 0, 0, 0.0007109993 ,0, 0, 0.003309997, 0.07398107, 0.08974432, 0.07314944, 0.0622156, 0.03105934, 0.010014, 0.001267644, 0 ,0 ,0 ,0 ,0.020851, 0, 0, 0, 0, 0, 0, 0,
# 0, 0 ,0 ,0 ,0 ,0 ,0.0008845589, 0.00388914, 0, 0, 0, 0, 0.004711973, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0 ,0 ,0 ,0,
# 0, 0, 0, 0, 0, 0 ,0.0006980912, 0.02864456, 0.0008348122, 0 ,0, 0 ,0.3265275 ,0.01410418 ,0 ,0 ,0, 0.332281, 0, 0.001184787, 0, 0 ,0 ,0 ,0,
# 0, 0, 0, 0, 0, 0, 0.001044587, 0.005987847, 0.002030869 ,0, 0 ,0 ,0.0419751, 0.01410418, 0 ,0 ,0 ,0.027689 ,0, 0.0006171872, 0 ,0, 0, 0, 0,
# 0, 0, 0 ,0 ,0, 0, 0.002290674, 0.09854721, 0.1947723, 0.001041604, 0.0008129025, 0.013273, 0.1380028, 0.02820811, 0.015272, 0 ,0 ,0, 0.070001, 0.0006171872, 0, 0, 0, 0, 0,
# 0, 0, 0, 0, 0, 0.01216899, 0.004603415, 0.09481322, 0.001261869, 0.003441011, 0, 0, 0.281883, 0.08955499, 0 ,0 ,0,0.304589, 0.070001, 0.03959263, 0, 0, 0, 0, 0,
# 0, 0, 0, 0, 0, 5.899995E-05, 0.08581191, 0.008719152, 0.001671519, 0.0001719213, 0.001211964, 0.001571, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
# 0.726929, 0, 0, 0.315373, 0, 0, 0.4453562, 0.5160776, 0.3443648, 0.8572081, 0.3833565, 0.890893, 0, 0 ,0.123605, 0 ,0 ,0 ,0 ,0.005774691, 0.161448, 0.055941, 0, 0 ,0,
# 0.050637, 0 ,0, 0, 0, 0, 0.01242108, 0.0372947, 0.09778056, 0.05176783, 0.5141716, 0, 0, 0, 0.244326, 0 ,0 ,0 ,0 ,0.005284418 ,0.696038 ,0.47728, 0.058572 ,0, 0,
# 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00386473, 0.095929, 0.11221, 0.161911, 0.099731, 0,
# 0, 0, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0 ,0 ,0 ,0.006854988, 0.008578, 0.04999995, 0.095, 0, 0, 0.0102941, 0.005032, 0.015172, 0.076645, 0.202492, 0,
# 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.002426699, 0, 0, 0, 0.069225, 0.4222216, 0.882656, 0 ,0, 0.08967555, 0.041553, 0.127274, 0.702872, 0.209335, 0,
# 0 ,0.005316995, 0, 0.0005529983, 0, 0, 0, 0, 0, 0, 0, 0, 0.002525342, 0.001344179, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0, 0, 0,
# 0 ,0 ,0, 0 ,0, 0, 0, 0, 0, 0 ,0, 0, 0.06347585, 0.8361391, 0.534461, 0.5277785, 0.022344, 0.2769, 0.839998, 0.8429105, 0, 0.212123, 0, 0.488442, 1,
# 0, 0, 0, 0.2077324, 0, 0.000843999, 0.009141359, 0.01302205, 0, 0, 0.006723996, 0, 0.03762027, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
# 29,25, byrow = T))
#----------------------------
# code to copy paste data into Rpath, can do this or copy paste to CSV and import.
# read.excel <- function(header=TRUE,...) {
# read.table("clipboard",sep="\t",header=F,...)
# }
#
# dcdat<-read.excel()
# DC.data<-as.data.table(dcdat)
# save(DC.data, file='DC.data.noscallop.RData')
load('DC.data.noscallop.RData')
#-------------------------------------------
#Fix column names
# setnames(WSS28.params$diet, 1, 'Group')
# setnames(WSS28.params$diet, paste0('V', 1:25), DC.groups)
# setnames(WSS28.params$diet, c(paste('V', 1:27, sep = '')), c('Group', DC.groups)
#Merge DC matrix with diet parameter object
WSS28.params$diet <- cbind(WSS28.params$diet[, Group], DC.data)
#Fix column names
setnames(WSS28.params$diet, 1, 'Group')
setnames(WSS28.params$diet, paste0('V', 1:25), DC.groups)
#Ecopath
WSS28 <- rpath(WSS28.params, 'Western Scotian Shelf')
save(WSS28.params, file = "WSS28.RData")
Make CModels WSS28
#QPress-------------------------------------------------------------------------
setwd("~/JCT work/WGNARS/RpathQNM/RpathQNM/WSS28model")
#Prey matrix
qpress.prey <- WSS28.params$diet
pred.cols <- names(qpress.prey)[which(names(qpress.prey) != 'Group')]
qpress.prey <- qpress.prey[Group != 'Import']
#Ensure that the matrices are square
detritus.DC <- data.table(Group = qpress.prey[, Group], Detritus = 0, Discards = 0)
qpress.prey <- merge(qpress.prey, detritus.DC, by = 'Group')
#line columns and rows up
groupnames <- as.character(qpress.prey$Group)
groupnames<-c('Group', groupnames)
setcolorder(qpress.prey, groupnames)
#Predator M2/F matrix
x <- copy(WSS28)
ngroup <- x$NUM_LIVING + x$NUM_DEAD
out <- data.frame(Group = x$Group[1:ngroup],
type = x$type [1:ngroup],
PB = x$PB [1:ngroup])
#Calculate F mortality
# totcatch <- x$Catch + x$Discards
totcatch <- x$Landings + x$Discards
Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))])
# setnames(Fmort, 'Fleet1', 'Fishery')
# Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))])
# setnames(Fmort, paste('V', 1:x$NUM_GEARS, sep = ''),
# paste('F.', x$Group[(ngroup +1):x$NUM_GROUPS], sep = ''))
out <- cbind(out, Fmort[1:ngroup, ])
setnames(out, 'Fmort[1:ngroup, ]', 'Fishery')
# pred.cols <- names(qpress.prey)[which(names(qpress.prey) != 'Group')]
#
# qpress.prey <- qpress.prey[Group != 'Import']
#
# #Ensure that the matrices are square
# detritus.DC <- data.table(Group = qpress.prey[, Group], Detritus = 0, Discards = 0)
# qpress.prey <- merge(qpress.prey, detritus.DC, by = 'Group')
#
# #line columns and rows up
# groupnames <- as.character(qpress.prey$Group)
# groupnames<-c('Group', groupnames)
# setcolorder(qpress.prey, groupnames)
#
# #Predator M2/F matrix
# x <- copy(WSS28)
# ngroup <- x$NUM_LIVING + x$NUM_DEADs
# out <- data.frame(Group = x$Group[1:ngroup],
# type = x$type [1:ngroup],
# PB = x$PB [1:ngroup])
# #Calculate F mortality
# totcatch <- x$Catch + x$Discards
# Fmort <- as.data.frame(totcatch / x$BB[row(as.matrix(totcatch))])
# setnames(Fmort, paste('V', 1:x$NUM_GEARS, sep = ''),
# paste('F.', x$Group[(ngroup +1):x$NUM_GROUPS], sep = ''))
# out <- cbind(out, Fmort[1:ngroup, ])
# setnames(out, 'Fmort[1:ngroup, ]', 'Fishery')
#Add Fishery to prey
fishery <- as.data.table(out)
fishery <- fishery[, list(Group, Fishery)]
qpress.prey <- merge(qpress.prey, fishery, by="Group")
qpress.prey<-qpress.prey[, -1]#remove Group column
fishery.0 <- as.data.table(matrix(rep(0, nrow(fishery) + 1), 1, nrow(fishery) + 1))
# setnames(fishery.0, paste('V', 1:x$NUM_Groups, sep=''), groupnames)
qpress.prey <- rbindlist(list(qpress.prey, fishery.0), use.names=FALSE)
# fishery <- as.data.table(out)
# fishery <- fishery[, list(Group, Fishery)]
# qpress.prey <- merge(qpress.prey, fishery, by = 'Group')
# qpress.prey[, Group := NULL]
# fishery.0 <- as.data.table(matrix(rep(0, nrow(fishery) + 1), 1, nrow(fishery) + 1))
# qpress.prey <- rbindlist(list(qpress.prey, fishery.0))
# #Predator M2/F matrix
# x <- copy(WSS28.b)
# ngroup <- x$NUM_LIVING + x$NUM_DEAD
# out <- data.frame(Group = x$Group[1:ngroup],
# type = x$type [1:ngroup],
# PB = x$PB [1:ngroup])
#Calculate M0
M0 <- c(x$PB[1:x$NUM_LIVING] * (1 - x$EE[1:x$NUM_LIVING]),
x$EE[(x$NUM_LIVING + 1):ngroup])
out <- cbind(out, M0)
#Calculate M2
bio <- x$Biomass[1:x$NUM_LIVING]
BQB <- bio * x$QB[1:x$NUM_LIVING]
diet <- as.data.frame(x$DC)
nodetrdiet <- diet[1:x$NUM_LIVING, ]
detrdiet <- diet[(x$NUM_LIVING +1):ngroup, ]
newcons <- nodetrdiet * BQB[col(as.matrix(nodetrdiet))]
predM <- newcons / bio[row(as.matrix(newcons))]
detcons <- detrdiet * BQB[col(as.matrix(detrdiet))]
predM <- rbind(predM, detcons)
# setnames(predM, paste('V', 1:x$NUM_LIVING, sep = ''),
# paste('M2.', x$Group[1:x$NUM_LIVING], sep = ''))
qpress.pred <- as.data.table(cbind(out, predM))
M2.cols <- names(qpress.pred)[which(!names(qpress.pred) %in% c('Group', 'type',
'PB', 'M0'))]
qpress.pred[, mort.sum := rowSums(.SD), .SDcols = M2.cols]
qpress.pred <- qpress.pred[, .SD / mort.sum, .SDcols = M2.cols, by = Group]
qpress.pred[is.na(qpress.pred)] <- 0
#Code for thresholds
# for(ipred in seq_along(M2.cols)){
# setnames(qpress.pred, M2.cols[ipred], 'V1')
# qpress.pred[V1 < .1, V1 := 0]
# setnames(qpress.pred, 'V1', M2.cols[ipred])
# }
#make qpress.pred negative
qpress.pred <- qpress.pred[, .SD * -1, .SDcol = M2.cols, by = Group]
#Ensure that the matrices are square
detritus.m2 <- data.table(qpress.pred[, 'Group'], M2.Detritus = 0, M2.Discards = 0)
qpress.pred <- merge(qpress.pred, detritus.m2, by = 'Group')
qpress.pred <- qpress.pred[, .SD, .SDcols = c(M2.cols, 'M2.Detritus', 'M2.Discards')]
prednames <- as.character(c('Fishery', detritus.DC$Group))
colnames(qpress.pred) <- prednames
new.names<-colnames(qpress.prey)
setcolorder(qpress.pred, new.names)
qpress.pred <- t(qpress.pred)
qpress.pred <- cbind(qpress.pred, rep(0, nrow(qpress.pred)))
colnames(qpress.pred) <- c(new.names)
#Merge pred and prey matrices
qpress.WSS28 <- qpress.pred + qpress.prey
#add/change row and column names
FG<-as.vector(new.names)
qpress.WSS28<-cbind(FG, qpress.WSS28)
save(qpress.WSS28, file="qpress.WSS28.RData")
#remove values below 10%, 20%, 30%, 40%, 50%
WSS28qnm10<-qpress.WSS28
WSS28qnm10[WSS28qnm10 >=-0.10 & WSS28qnm10<= 0.10]<-0
WSS28qnm20<-qpress.WSS28
WSS28qnm20[WSS28qnm20 >=-0.20 & WSS28qnm20<= 0.20]<-0
WSS28qnm30<-qpress.WSS28
WSS28qnm30[WSS28qnm20 >=-0.30 & WSS28qnm20<= 0.30]<-0
WSS28qnm40<-qpress.WSS28
WSS28qnm40[WSS28qnm40 >=-0.40 & WSS28qnm40<= 0.40]<-0
WSS28qnm50<-qpress.WSS28
WSS28qnm50[WSS28qnm50 >=-0.50 & WSS28qnm20<= 0.50]<-0
WSS28qnm100<-qpress.WSS28
save(WSS28qnm10, file='WSS28qnm10.RData')
save(WSS28qnm20, file='WSS28qnm20.RData')
save(WSS28qnm30, file='WSS28qnm30.RData')
save(WSS28qnm40, file='WSS28qnm40.RData')
save(WSS28qnm50, file='WSS28qnm50.RData')
save(WSS28qnm100, file='WSS28qnm100.RData')
Adjacency to digraph functions
MM2Qpress <- function(data){
library(data.table)
mental.sheet <- as.data.table(data)
names(mental.sheet)[1] <- 'Box'
n <- nrow(mental.sheet)
box.names <- names(mental.sheet)[which(names(mental.sheet) != 'Box')]
model <- c()
for(i in 1:n){
pos <- which(mental.sheet[i, .SD, .SDcol = box.names] > 0)
if(length(pos) > 0){
pos.interaction <- paste(mental.sheet[i, Box], '->',
names(mental.sheet[i, pos + 1, with = F]),
sep = '')
model <- append(model, pos.interaction)
}
neg <- which(mental.sheet[i, .SD, .SDcol = box.names] < 0)
if(length(neg) > 0){
neg.interaction <- paste(mental.sheet[i, Box], '-*',
names(mental.sheet[i, neg + 1, with = F]),
sep = '')
model <- append(model, neg.interaction)
}
}
return(model)
}
#function to create qnm models for qpress from signed digraphs with self-limitation
make.qnm<-function(modmat){
q<-MM2Qpress(modmat)
qnm<-dput(q)
qnm.mod<-parse.digraph(qnm)
qnm.model<-enforce.limitation(qnm.mod)
}
#function to create qnm models for qpress from signed digraphs WITHOUT self-limitation
make.qnm.nolim<-function(modmat){
q<-MM2Qpress(modmat)
qnm<-dput(q)
qnm.mod<-parse.digraph(qnm)
qnm.model<-qnm.mod
}
Qpress for WSS28
#--------------------------------------------------------------------
#####PERTURBATIONS#####
#to loop over multiple models
#create folder with files in directory, call list of files
set.seed(1984)
#load RData files if required
#create list
mod<-list(WSS28qnm10, WSS28qnm20, WSS28qnm30, WSS28qnm40, WSS28qnm50, WSS28qnm100)
#get MM matrix into qpress format, create model
q.models<-vector('list', length(mod))
for(i in seq_along(mod)){
q.models[[i]]<-make.qnm(mod[[i]])
}
#write dia files
# for(i in seq_along(q.models)){
# #can add text below to further specify name
# write.dia(q.models[[i]], paste0("WSS", i,".dia"))
# }
#loop through models to create results csv file and pdf barplots
for(j in seq_along(q.models)){
A <- adjacency.matrix(q.models[[j]])
#A
## Function to generate the community matrix
s <- community.sampler(q.models[[j]])
## Function to check the validation condition
#press <- press.validate(GB2el,
# perturb=c("Forage_Fish"=-1),
# monitor=c("Commercial_Pelagic_Fishery"=1))
#Press/Perturbation with decrease forage fish and no monitoring)
press <- press.validate(q.models[[j]],
perturb=c("Seals"=1), monitor=F)
## Function to define the perturbation scenario
impact <- press.impact(q.models[[j]],perturb=c("Seals"=1))
## Use 1000 simulations
n.sims <- 1000
results <- 0
i <- 0
while(i < n.sims) {
## Randomly choose edges to retain
z <- s$select(runif(1))
## Sample community matrix
W <- s$community()
## Check press condition and stability
if(!(press(W) && stable.community(W))) next
## Monitor impact post press
imp <- impact(W)
imp[abs(imp)<1e-07]=0 #added to remove small perturbations
results <- results + outer(sign(imp),-1:1,'==')
i <- i+1
rownames(results) <- levels(q.models[[j]]$From)
colnames(results) <- c('-','0','+')
save(results, file=paste0("WSS28","results", j,"Seals_up",".RData"))
write.csv(results, paste0("WSS28","results", j,"Seals_up",".csv"))
# library(RColorBrewer)
# pal <- brewer.pal(n=5,"Greys")[5:1]
# #opar <- par(mar=c(5,10,1,1)+0.1)
# prop <- results/rowSums(results)
# r <- colSums(t(prop)*(-1:1))
# filename= paste("WSS28", "plots", j,"Seals_up", ".pdf")
# pdf(filename)
# #adjust barplot if you want it ordered by proportion
# par(mar=c(5,10,1,1)+0.1)
# barplot(t(prop),
# horiz=T,cex.names=0.7,cex.axis=0.8,las=2,border=T, col=pal,xlab="Proportion")
if(dev.cur() > 1) dev.off()
# dev.off()
}
}
#code for ordered barplot
# pal <- brewer.pal(n=5,"Greys")[5:1]
# opar <- par(mar=c(5,10,1,1)+0.1)
# prop <- results/rowSums(results)
# r <- colSums(t(prop)*(-1:1))
# filename= paste("wss", "plot", j,"phytoplankton up", ".pdf")
# pdf(filename)
# barplot(t(prop[order(r),]),
# horiz=T,cex.names=0.6,cex.axis=0.8,las=2,border=T,col=pal,xlab="Proportion")
# par(opar)
Qpress for WSS28 no self limiting
#--------------------------------------------------------------------
#####PERTURBATIONS#####
#to loop over multiple models
#create folder with files in directory, call list of files
set.seed(1984)
#load RData files if required
#create list
mod<-list(WSS28qnm10, WSS28qnm20, WSS28qnm30, WSS28qnm40, WSS28qnm50, WSS28qnm100)
#get MM matrix into qpress format, create model
q.models<-vector('list', length(mod))
for(i in seq_along(mod)){
q.models[[i]]<-make.qnm.nolim(mod[[i]])
}
#write dia files
# for(i in seq_along(q.models)){
# #can add text below to further specify name
# write.dia(q.models[[i]], paste0("WSS", i,".dia"))
# }
#loop through models to create results csv file and pdf barplots
for(j in seq_along(q.models)){
A <- adjacency.matrix(q.models[[j]])
#A
## Function to generate the community matrix
s <- community.sampler(q.models[[j]])
## Function to check the validation condition
#press <- press.validate(GB2el,
# perturb=c("Forage_Fish"=-1),
# monitor=c("Commercial_Pelagic_Fishery"=1))
#Press/Perturbation with decrease forage fish and no monitoring)
press <- press.validate(q.models[[j]],
perturb=c("Phytoplankton"=1), monitor=F)
## Function to define the perturbation scenario
impact <- press.impact(q.models[[j]],perturb=c("Phytoplankton"=1))
## Use 1000 simulations
n.sims <- 1000
results <- 0
i <- 0
while(i < n.sims) {
## Randomly choose edges to retain
z <- s$select(runif(1))
## Sample community matrix
W <- s$community()
## Check press condition and stability
if(!(press(W) && stable.community(W))) next
## Monitor impact post press
imp <- impact(W)
imp[abs(imp)<1e-07]=0 #added to remove small perturbations
results <- results + outer(sign(imp),-1:1,'==')
i <- i+1
rownames(results) <- levels(q.models[[j]]$From)
colnames(results) <- c('-','0','+')
save(results, file=paste0("WSS28nolim","results", j,"Seals_up",".RData"))
write.csv(results, paste0("WSS28nolim","results", j,"Seals_up",".csv"))
# library(RColorBrewer)
# pal <- brewer.pal(n=5,"Greys")[5:1]
# #opar <- par(mar=c(5,10,1,1)+0.1)
# prop <- results/rowSums(results)
# r <- colSums(t(prop)*(-1:1))
# filename= paste("WSS28", "plots", j,"Seals_up", ".pdf")
# pdf(filename)
# #adjust barplot if you want it ordered by proportion
# par(mar=c(5,10,1,1)+0.1)
# barplot(t(prop),
# horiz=T,cex.names=0.7,cex.axis=0.8,las=2,border=T, col=pal,xlab="Proportion")
if(dev.cur() > 1) dev.off()
# dev.off()
}
}
#code for ordered barplot
# pal <- brewer.pal(n=5,"Greys")[5:1]
# opar <- par(mar=c(5,10,1,1)+0.1)
# prop <- results/rowSums(results)
# r <- colSums(t(prop)*(-1:1))
# filename= paste("wss", "plot", j,"phytoplankton up", ".pdf")
# pdf(filename)
# barplot(t(prop[order(r),]),
# horiz=T,cex.names=0.6,cex.axis=0.8,las=2,border=T,col=pal,xlab="Proportion")
# par(opar)
visualize results
WSS28_10_phyto_up<-get(load("WSS28results1Phytoplankton_up.RData"))
WSS28_10_phyto_down<-get(load("WSS28results1Phytoplankton_dtesting.RData"))
WSS28_20_phyto_up<-get(load("WSS28results2Phytoplankton_up.RData"))
WSS28_20_phyto_down<-get(load("WSS28results2Phytoplankton_dtesting.RData"))
WSS28_30_phyto_up<-get(load("WSS28results3Phytoplankton_up.RData"))
WSS28_30_phyto_down<-get(load("WSS28results3Phytoplankton_dtesting.RData"))
WSS28_40_phyto_up<-get(load("WSS28results4Phytoplankton_up.RData"))
WSS28_40_phyto_down<-get(load("WSS28results4Phytoplankton_dtesting.RData"))
WSS28_50_phyto_up<-get(load("WSS28results5Phytoplankton_up.RData"))
WSS28_50_phyto_down<-get(load("WSS28results5Phytoplankton_dtesting.RData"))
WSS28_100_phyto_up<-get(load("WSS28results6Phytoplankton_up.RData"))
WSS28_100_phyto_down<-get(load("WSS28results6Phytoplankton_dtesting.RData"))
WSS28_10_spelagics_up<-get(load("WSS28results1Small pelagics_up.RData"))
WSS28_10_spelagics_down<-get(load("WSS28results1Small pelagics_down.RData"))
WSS28_20_spelagics_up<-get(load("WSS28results2Small pelagics_up.RData"))
WSS28_20_spelagics_down<-get(load("WSS28results2Small pelagics_down.RData"))
WSS28_30_spelagics_up<-get(load("WSS28results3Small pelagics_up.RData"))
WSS28_30_spelagics_down<-get(load("WSS28results3Small pelagics_down.RData"))
WSS28_40_spelagics_up<-get(load("WSS28results4Small pelagics_up.RData"))
WSS28_40_spelagics_down<-get(load("WSS28results4Small pelagics_down.RData"))
WSS28_50_spelagics_up<-get(load("WSS28results5Small pelagics_up.RData"))
WSS28_50_spelagics_down<-get(load("WSS28results5Small pelagics_down.RData"))
WSS28_100_spelagics_up<-get(load("WSS28results6Small pelagics_up.RData"))
WSS28_100_spelagics_down<-get(load("WSS28results6Small pelagics_down.RData"))
WSS28_10_Seals_up<-get(load("WSS28results1Seals_up.RData"))
WSS28_10_Seals_down<-get(load("WSS28results1Seals_down.RData"))
WSS28_20_Seals_up<-get(load("WSS28results2Seals_up.RData"))
WSS28_20_Seals_down<-get(load("WSS28results2Seals_down.RData"))
WSS28_30_Seals_up<-get(load("WSS28results3Seals_up.RData"))
WSS28_30_Seals_down<-get(load("WSS28results3Seals_down.RData"))
WSS28_40_Seals_up<-get(load("WSS28results4Seals_up.RData"))
WSS28_40_Seals_down<-get(load("WSS28results4Seals_down.RData"))
WSS28_50_Seals_up<-get(load("WSS28results5Seals_up.RData"))
WSS28_50_Seals_down<-get(load("WSS28results5Seals_down.RData"))
WSS28_100_Seals_up<-get(load("WSS28results6Seals_up.RData"))
WSS28_100_Seals_down<-get(load("WSS28results6Seals_down.RData"))
WSS28_10_Fishery_up<-get(load("WSS28results1Fishery_up.RData"))
WSS28_10_Fishery_down<-get(load("WSS28results1Fishery_down.RData"))
WSS28_20_Fishery_up<-get(load("WSS28results2Fishery_up.RData"))
WSS28_20_Fishery_down<-get(load("WSS28results2Fishery_down.RData"))
WSS28_30_Fishery_up<-get(load("WSS28results3Fishery_up.RData"))
WSS28_30_Fishery_down<-get(load("WSS28results3Fishery_down.RData"))
WSS28_40_Fishery_up<-get(load("WSS28results4Fishery_up.RData"))
WSS28_40_Fishery_down<-get(load("WSS28results4Fishery_down.RData"))
WSS28_50_Fishery_up<-get(load("WSS28results5Fishery_up.RData"))
WSS28_50_Fishery_down<-get(load("WSS28results5Fishery_down.RData"))
WSS28_100_Fishery_up<-get(load("WSS28results6Fishery_up.RData"))
WSS28_100_Fishery_down<-get(load("WSS28results6Fishery_down.RData"))
WSS28.qnm.list<-list()
files <- list.files(full.names = T, recursive = T, pattern = "biomass_annual.csv")
subfiles<-files[c(T,F)]
nfiles<-length(subfiles)
#create empty dataframe with rows 1:49 (change if using a different model with differnt number of FGs)
biomass_mc<-data.frame()[1:49,]
for(i in 1:nfiles){
biomass<-read.csv(paste(subfiles[[i]]),skip=9,header=TRUE)
#pulls data from year 1, change by adjusting row number
bio_y<-biomass[2,-1]
rownames(bio_y)<-(paste0("X",i))
biomass_t<-t(bio_y)
biomass_mc<-cbind(biomass_mc, biomass_t)
}
single Qpress for WSS28
data<-WSS28qnm50
WSS28qnm50<-MM2Qpress(data)
WSS28qnm50<-dput(WSS28qnm50)
WSS28qnm50<-parse.digraph(WSS28qnm50)
# write.dia(WSS28qnm10Model, file="WSS28test.dia", width=8, height=2, self=T)
#enforce
# WSS28qnm100el<-enforce.limitation(WSS28qnm100)
## Examine unweighted adjacency matrix
A <- adjacency.matrix(WSS28qnm50)
# A <- adjacency.matrix(WSS28qnm50el)
## Function to generate the community matrix
s <- community.sampler(WSS28qnm50)
# s <- community.sampler(WSS28qnm100el)
## Function to check the validation condition
#press <- press.validate(GB2el,
# perturb=c("Forage_Fish"=-1),
# monitor=c("Commercial_Pelagic_Fishery"=1))
#Press/Perturbation with decrease forage fish and no monitoring)
press <- press.validate(WSS28qnm50,
perturb=c("Phytoplankton"=1), monitor=F)
## Function to define the perturbation scenario
impact <- press.impact(WSS28qnm50,perturb=c("Phytoplankton"=1))
## Use 10000 simulations
n.sims <- 1000
results <- 0
i <- 0
while(i < n.sims) {
## Randomly choose edges to retain
z <- s$select(runif(1))
## Sample community matrix
W <- s$community()
## Check press condition and stability
if(!(press(W) && stable.community(W))) next
## Monitor impact post press
imp <- impact(W)
imp[abs(imp)<1e-07]=0 #added to remove small perturbations
results <- results + outer(sign(imp),-1:1,'==')
i <- i+1
}
## Print results
rownames(results) <- levels(WSS28qnm50$From)
colnames(results) <- c('-','0','+')
#results
##save results
WSS28qnm50_phyto_down<-results
write.csv(WSS28qnm50_phyto_down, file="WSS2850nolimit_phyto_up.csv")
save(WSS28qnm50_phyto_down, file="WSS2850nolimit_phyto_up.RData")
Rpath for WSS62
## Western Scotian Shelf model
#-------------------------------------------------------------------------------
#Required packages
library(data.table); library(Rpath); library(here)
#-------------------------------------------------------------------------------
#User created functions
#-------------------------------------------------------------------------------
#Western Scotian Shelf Groups
groups <- c('Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks',
'Large pelagic', 'Cod <1', 'Cod 1-3', 'Cod 4-6', 'Cod 7+',
'S Hake <25', 'S Hake 25-31', 'S Hake 31+', 'Halibut <46', 'Halibut 46-81',
'Halibut 82+', 'Pollock <49', 'Pollock 49+', 'D piscvores <40', 'D piscvores 40+',
'L benthivores <41', 'L bentivores 41+', 'Skates <49', 'Skates 49+',
'Dogfish', 'Redfish <22', 'Redfish 22+', 'A plaice <26', 'A plaice 26+',
'Flounders <30', 'Flounders 30+', 'Haddock <3', 'Haddock 3+', 'L sculpin <25',
'L sculpin 25+', 'Herring <4', 'Herring 4+', 'Other pelagic', 'Mackerel',
'Mesopelagic', 'Small-medium benthivores', 'Squids', 'Lobster', 'Large crabs',
'Small crabs', 'Shrimps', 'Scallop', 'Bivalves', 'Other molluscs',
'Other arthropoda', 'Echinoderms', 'Sessile benthic groups', 'Worms',
'Meiofauna','Gelatinous zoop', 'Macrozoop', 'Mesozoop', 'Microzoop',
'Microflora', 'Phytoplankton', 'Discards', 'Detritus', 'Fleet1')
#Identify type (0 = consumer, 1 = producer, 2 = detritus, 3 = fleet)
types <- c(rep(0, 59), 1, 2, 2, 3)
#Identify stanza groups
stgroups <- c(rep(NA, 6), rep('Cod', 4), rep('S hake', 3), rep('Halibut', 3),
rep('Pollock', 2), rep('D piscvores', 2), rep('L benthivores', 2),
rep('Skates', 2), NA, rep('Redfish', 2), rep('A plaice', 2),
rep('Flounders', 2), rep('Haddock', 2), rep('S sculpin', 2),
rep('Herring', 2), rep(NA, 26))
#Create parameter list object
WSS62.params <- create.rpath.params(groups, types, stgroups)
#Fill in model parameters
#Biomass, Production, consumption
biomass <- c(0.407075, 0.049754, 0.04383, 0.00617, 0.0260249, 0.023573, rep(NA, 3),
0.076371, rep(NA, 2), 0.0473429, rep(NA, 2), 0.017291, NA, 0.154529,
NA, 0.26641, NA, 0.071238, NA, 0.073658, 2.102007, NA, 1.523083, NA,
0.101483, NA, 0.25631, NA, 0.806943, NA, 0.073649, NA, 2.25021, 0.86577,
0.535679, 0.012515, 0.037782, 0.173893, 0.301295, 0.157303, 0.965481,
1.302563, 0.646857, 63.37292, 2, 0.853, 12.63734, 26.02519, 7.345485,
4.09129, 0.5204, 41.31, 23.79555, 5.800508, 3.678861, 33.664, 0.063293,
1, NA)
pb <- c(0.071, 0.18, 0.147007, 0.25, 0.18, 0.4, 1.55, 0.39, 0.82, 0.82, 1.6, 1.3,
1.3, 0.57, 0.27, 0.27, 0.47, 0.85, 0.44, 0.62, 0.62, 0.64, 0.2, 0.26, 0.1393199,
0.27, 0.23, 0.44, 0.52, 0.53, 0.59, 0.66, 0.42, 0.6, 0.5, 0.64, 0.81, 0.74,
0.583, 0.97, 1.41, 4, 0.91, 0.654, 1.31, 3, 0.342, 0.69, 0.75, 4.4, 0.33,
0.32, 1.25, 10.8, 15.51, 3.04, NA, 82.8, 104.938, 70.639, rep(NA, 3))
qb <- c(4.94, 14.5, 7.338632, 87.6, 4.78, 4.24, rep(NA, 3), 1.326, rep(NA, 2), 3.28,
rep(NA, 2), 1.61, NA, 2.94, NA, 2.5, NA, 1.92, NA, 1.9, 1.739802, NA, 2.41,
NA, 2.41, NA, 3.21, NA, 2.08, NA, 3.93, NA, 2.36, 3.31, 5.08, 18.9, 4.7,
11.33333, rep(NA, 12), 62.05, 19.5, 73, rep(NA, 6))
WSS62.params$model[, Biomass := biomass]
WSS62.params$model[, PB := pb]
WSS62.params$model[, QB := qb]
#Production to Consumption
WSS62.params$model[Group %in% c('Lobster', 'Large crabs', 'Small crabs', 'Shrimps',
'Scallop', 'Bivalves', 'Other molluscs', 'Other arthropoda',
'Echinoderms', 'Sessile benthic groups', 'Worms',
'Meiofauna'),
ProdCons := 0.15]
WSS62.params$model[Group == 'Mesozoop', ProdCons := 0.4]
WSS62.params$model[Group %in% c('Microzoop', 'Microflora'), ProdCons := 0.5]
#Biomass accumulation and unassimilated production
# ba <- c(0.003257, 0, 0.003702, 0, -0.001952, 0, -0.000180536, -0.0111684, -0.01359565,
# -0.00381855, rep(0, 6), -0.005120684, -0.00618116, -0.00398017, -0.0133205,
# -0.008473957, -0.0071238, rep(0, 5), -0.002998298, -0.00710381, 0, 0,
# 0.02020396, 0.05648601, -0.003061481, -0.00368245, rep(0, 7), 0.0301295,
# rep(0, 3), 0.03234285, rep(0, 13), rep(NA, 3))
#Needed to turn off biomass accumulation for stanza groups as it is already baked in
ba <- c(0.003257, 0, 0.003702, 0, -0.001952, rep(0, 37), 0.0301295,
rep(0, 3), 0.03234285, rep(0, 13), rep(NA, 3))
WSS62.params$model[, BioAcc := ba]
WSS62.params$model[, Unassim := c(rep(0.2, 55), rep(0.4, 3), 0.2, rep(0, 3), NA)]
#Detrital fate
WSS62.params$model[, Detritus := c(rep(1, 60), rep(0, 3))]
WSS62.params$model[, Discards := c(rep(0, 62), 1)]
#Landings/Discards
land <- c(rep(0, 4), 0.00209077, 0.005067847, 0, 0.024145, 0.08195, 0.010009,
0.005438, 0.023001, 0.008442, 2.9E-05, 0.00025, 0.00299, 0.004407,
0.092313, 0, 0.058946, 0, 0.004406, 0, 0.000799, 0.01325133, 0,
0.05914, 0, 0.005274, 0, 0.021106, 0.000433, 0.089896, 0.00084, 0.001181,
0.230954, 0.719899, 0.019708, 0.048208, 0, 0, 0.001637, 0.203756,
0.016668, 0, 0.000162, 0.106389, 0.010353, 0.00236, 0, 0.03026,
rep(0, 11), NA)
disc <- c(0, 0.002995, 0, 2.4E-05, 0.000575, 0.000469, 0.000286, rep(0, 3), 9.6E-05,
0, 0, 3.2E-05, 0.000274, 0, 0.00029, 0, 0.002805, 0, 0.000385, 0,
0.003928, 0.009166, 0.009694, 0.00027, 0, 9E-06, 0, 0.000372, 0, 0.000111,
0, 0.000367, 0.000516, 0.000211, 0, 0.00056, 0.000276, 0, 0.000236,
1.9E-05, 0.001526, 0.00725, 0.000633, 8E-06, 0, 0, 4.1E-05, 0, 0.016926,
0.002943, rep(0, 10), NA)
WSS62.params$model[, Fleet1 := land]
WSS62.params$model[, Fleet1.disc := disc]
#Diet composition
#read DC from EwE6.6
DC.groups <- groups[which(!groups %in% c('Discards', 'Detritus', 'Fleet1'))]
DC.data <- as.data.table(read.csv(here('data-raw', 'WSS and BOF 11-Diet composition.csv')))
setnames(DC.data, names(DC.data), c('X', 'Group', DC.groups))
DC.data[, c('X', 'Phytoplankton') := NULL]
DC.data <- DC.data[1:63, ]
DC.data[, Phytoplankton := as.numeric(0)]
#Merge DC matrix with diet parameter object
# WSS.params$diet <- cbind(WSS.params$diet[, Group], DC.data)
# #Fix column names
# setnames(WSS.params$diet, 1, 'Group')
# setnames(WSS.params$diet, paste0('V', 1:60), DC.groups)
WSS62.params$diet <- DC.data
#Set multistanza parameters
#Group parameters
WSS62.params$stanzas$stgroups[, VBGF_Ksp := c(0.14, 0.46, 0.14, 0.16, 0.075, 0.14,
0.13, 0.06, 0.12, 0.21, 0.27, 0.12,
0.28)]
WSS62.params$stanzas$stgroups[, Wmat := c(0.04, 0.22, 0.24, 0.13, 0.02, 0.15,
0.4, 0.06, 0.06, 0.17, 0.16, 0.15,
0.33)]
WSS62.params$stanzas$stgroups[, BAB := c(-0.05, 0, 0, -0.04, -0.05, -0.10, 0,
0, -0.07, 0, 0.07, -0.05, 0)]
#Individual stanza parameters
WSS62.params$stanzas$stindiv[, First := c(0, 12, 48, 84, 0, 25, 48, 0, 36, 72,
0, 48, 0, 47, 0, 66, 0, 109, 0, 108,
0, 57, 0, 46, 0, 36, 0, 65, 0, 48)]
WSS62.params$stanzas$stindiv[, Last := c(11, 47, 83, 400, 24, 47, 400, 35, 71,
400, 47, 400, 46, 400, 65, 400, 108,
400, 107, 400, 56, 400, 45, 400, 35,
400, 64, 400, 47, 400)]
WSS62.params$stanzas$stindiv[, Z := c(1.55, 0.39, 0.82, 0.82, 1.6, 1.3, 1.3,
0.57, 0.27, 0.27, 0.47, 0.85, 0.44, 0.62,
0.62, 0.64, 0.2, 0.26, 0.27, 0.23, 0.44,
0.52, 0.53, 0.59, 0.66, 0.42, 0.6, 0.5,
0.64, 0.81)]
WSS62.params$stanzas$stindiv[, Leading := c(F, F, F, T, F, F, T, F, F, T, rep(c(F, T), 10))]
#Need to run this function to calculate all the stanza biomass/consumptions
#The calculated biomasses are different than EwE because of the BA terms...will investigate further
#Removed BioAcc for stanza groups and this works close enough (rounding errors)
WSS62.params <- rpath.stanzas(WSS62.params)
#Ecopath
WSS62 <- rpath(WSS62.params, 'Western Scotian Shelf')
usethis::use_data(WSS62, overwrite = TRUE)
Make CModels WSS62
#QPress-------------------------------------------------------------------------
setwd("~/JCT work/WGNARS/RpathQNM/RpathQNM/WSS62model")
Group <- c( 'Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks',
'Large pelagic', 'Cod less1', 'Cod 1to3', 'Cod 4to6', 'Cod 7plus',
'S Hake less25', 'S Hake 25to31', 'S Hake 31plus', 'Halibut less46', 'Halibut 46to81',
'Halibut 82plus', 'Pollock less49', 'Pollock 49plus', 'D piscvores less40', 'D piscvores 40plus',
'L benthivores less41', 'L bentivores 41plus', 'Skates less49', 'Skates 49plus',
'Dogfish', 'Redfish less22', 'Redfish 22plus', 'A plaice less26', 'A plaice 26plus',
'Flounders less30', 'Flounders 30plus', 'Haddock less3', 'Haddock 3plus', 'L sculpin less25',
'L sculpin 25plus', 'Herring less4', 'Herring 4plus', 'Other pelagic', 'Mackerel',
'Mesopelagic', 'Smallmedium benthivores', 'Squids', 'Lobster', 'Large crabs',
'Small crabs', 'Shrimps', 'Scallop', 'Bivalves', 'Other molluscs',
'Other arthropoda', 'Echinoderms', 'Sessile benthic groups', 'Worms',
'Meiofauna','Gelatinous zoop', 'Macrozoop', 'Mesozoop', 'Microzoop',
'Microflora', 'Phytoplankton', 'Discards', 'Detritus', 'Import')
groupnamescol<-c('Group', 'Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks',
'Large pelagic', 'Cod less1', 'Cod 1to3', 'Cod 4to6', 'Cod 7plus',
'S Hake less25', 'S Hake 25to31', 'S Hake 31plus', 'Halibut less46', 'Halibut 46to81',
'Halibut 82plus', 'Pollock less49', 'Pollock 49plus', 'D piscvores less40', 'D piscvores 40plus',
'L benthivores less41', 'L bentivores 41plus', 'Skates less49', 'Skates 49plus',
'Dogfish', 'Redfish less22', 'Redfish 22plus', 'A plaice less26', 'A plaice 26plus',
'Flounders less30', 'Flounders 30plus', 'Haddock less3', 'Haddock 3plus', 'L sculpin less25',
'L sculpin 25plus', 'Herring less4', 'Herring 4plus', 'Other pelagic', 'Mackerel',
'Mesopelagic', 'Smallmedium benthivores', 'Squids', 'Lobster', 'Large crabs',
'Small crabs', 'Shrimps', 'Scallop', 'Bivalves', 'Other molluscs',
'Other arthropoda', 'Echinoderms', 'Sessile benthic groups', 'Worms',
'Meiofauna','Gelatinous zoop', 'Macrozoop', 'Mesozoop', 'Microzoop',
'Microflora', 'Phytoplankton')
#Prey matrix
qpress.prey <- WSS62.params$diet
#change names to be compatible with Qpress
setnames(qpress.prey, groupnamescol)
qpress.prey<-qpress.prey[,-1]
qpress.prey<-cbind(qpress.prey, Group)
pred.cols <- names(qpress.prey)[which(names(qpress.prey) != 'Group')]
qpress.prey <- qpress.prey[Group != 'Import']
#Ensure that the matrices are square
detritus.DC <- data.table(Group = qpress.prey[, Group], Detritus = 0, Discards = 0)
qpress.prey <- merge(qpress.prey, detritus.DC, by = 'Group')
#line columns and rows up
groupnames <- as.character(qpress.prey$Group)
groupnames<-c('Group', groupnames)
setcolorder(qpress.prey, groupnames)
#Predator M2/F matrix
x <- copy(WSS62)
ngroup <- x$NUM_LIVING + x$NUM_DEAD
out <- data.frame(Group = x$Group[1:ngroup],
type = x$type [1:ngroup],
PB = x$PB [1:ngroup])
#Calculate F mortality
# totcatch <- x$Catch + x$Discards
totcatch <- x$Landings + x$Discards
Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))])
# setnames(Fmort, 'Fleet1', 'Fishery')
# Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))])
# setnames(Fmort, paste('V', 1:x$NUM_GEARS, sep = ''),
# paste('F.', x$Group[(ngroup +1):x$NUM_GROUPS], sep = ''))
out <- cbind(out, Fmort[1:ngroup, ])
setnames(out, 'Fmort[1:ngroup, ]', 'Fishery')
#Add Fishery to prey
fishery <- as.data.table(out)
fishery <- fishery[, list(Group, Fishery)]
qpress.prey <- merge(qpress.prey, fishery, by="Group")
qpress.prey<-qpress.prey[, -1]#remove Group column
fishery.0 <- as.data.table(matrix(rep(0, nrow(fishery) + 1), 1, nrow(fishery) + 1))
# setnames(fishery.0, paste('V', 1:x$NUM_Groups, sep=''), groupnames)
qpress.prey <- rbindlist(list(qpress.prey, fishery.0), use.names=FALSE)
# #Predator M2/F matrix
# x <- copy(WSS62.b)
# ngroup <- x$NUM_LIVING + x$NUM_DEAD
# out <- data.frame(Group = x$Group[1:ngroup],
# type = x$type [1:ngroup],
# PB = x$PB [1:ngroup])
#Calculate M0
M0 <- c(x$PB[1:x$NUM_LIVING] * (1 - x$EE[1:x$NUM_LIVING]),
x$EE[(x$NUM_LIVING + 1):ngroup])
out <- cbind(out, M0)
#Calculate M2
bio <- x$Biomass[1:x$NUM_LIVING]
BQB <- bio * x$QB[1:x$NUM_LIVING]
diet <- as.data.frame(x$DC)
nodetrdiet <- diet[1:x$NUM_LIVING, ]
detrdiet <- diet[(x$NUM_LIVING +1):ngroup, ]
newcons <- nodetrdiet * BQB[col(as.matrix(nodetrdiet))]
predM <- newcons / bio[row(as.matrix(newcons))]
detcons <- detrdiet * BQB[col(as.matrix(detrdiet))]
predM <- rbind(predM, detcons)
# setnames(predM, paste('V', 1:x$NUM_LIVING, sep = ''),
# paste('M2.', x$Group[1:x$NUM_LIVING], sep = ''))
qpress.pred <- as.data.table(cbind(out, predM))
M2.cols <- names(qpress.pred)[which(!names(qpress.pred) %in% c('Group', 'type',
'PB', 'M0'))]
qpress.pred[, mort.sum := rowSums(.SD), .SDcols = M2.cols]
qpress.pred <- qpress.pred[, .SD / mort.sum, .SDcols = M2.cols, by = Group]
qpress.pred[is.na(qpress.pred)] <- 0
#Code for thresholds
# for(ipred in seq_along(M2.cols)){
# setnames(qpress.pred, M2.cols[ipred], 'V1')
# qpress.pred[V1 < .1, V1 := 0]
# setnames(qpress.pred, 'V1', M2.cols[ipred])
# }
#make qpress.pred negative
qpress.pred <- qpress.pred[, .SD * -1, .SDcol = M2.cols, by = Group]
#Ensure that the matrices are square
detritus.m2 <- data.table(qpress.pred[, 'Group'], M2.Detritus = 0, M2.Discards = 0)
qpress.pred <- merge(qpress.pred, detritus.m2, by = 'Group')
qpress.pred <- qpress.pred[, .SD, .SDcols = c(M2.cols, 'M2.Detritus', 'M2.Discards')]
prednames <- as.character(c('Fishery', detritus.DC$Group))
colnames(qpress.pred) <- prednames
new.names<-colnames(qpress.prey)
setcolorder(qpress.pred, new.names)
qpress.pred <- t(qpress.pred)
qpress.pred <- cbind(qpress.pred, rep(0, nrow(qpress.pred)))
colnames(qpress.pred) <- c(new.names)
#Merge pred and prey matrices
qpress.WSS62 <- qpress.pred + qpress.prey
#add/change row and column names
FG<-as.vector(new.names)
qpress.WSS62<-cbind(FG, qpress.WSS62)
save(qpress.WSS62, file="qpress.WSS62.RData")
#remove values below 10%, 20%, 30%, 40%, 50%
WSS62qnm10<-qpress.WSS62
WSS62qnm10[WSS62qnm10 >=-0.10 & WSS62qnm10<= 0.10]<-0
WSS62qnm20<-qpress.WSS62
WSS62qnm20[WSS62qnm20 >=-0.20 & WSS62qnm20<= 0.20]<-0
WSS62qnm30<-qpress.WSS62
WSS62qnm30[WSS62qnm20 >=-0.30 & WSS62qnm20<= 0.30]<-0
WSS62qnm40<-qpress.WSS62
WSS62qnm40[WSS62qnm40 >=-0.40 & WSS62qnm40<= 0.40]<-0
WSS62qnm50<-qpress.WSS62
WSS62qnm50[WSS62qnm50 >=-0.50 & WSS62qnm20<= 0.50]<-0
WSS62qnm100<-qpress.WSS62
save(WSS62qnm10, file='WSS62qnm10.RData')
save(WSS62qnm20, file='WSS62qnm20.RData')
save(WSS62qnm30, file='WSS62qnm30.RData')
save(WSS62qnm40, file='WSS62qnm40.RData')
save(WSS62qnm50, file='WSS62qnm50.RData')
save(WSS62qnm100, file='WSS62qnm100.RData')
Qpress for WSS62
setwd("~/JCT work/WGNARS/RpathQNM/RpathQNM/WSS62model")
MM2Qpress <- function(data){
library(data.table)
mental.sheet <- as.data.table(data)
names(mental.sheet)[1] <- 'Box'
n <- nrow(mental.sheet)
box.names <- names(mental.sheet)[which(names(mental.sheet) != 'Box')]
model <- c()
for(i in 1:n){
pos <- which(mental.sheet[i, .SD, .SDcol = box.names] > 0)
if(length(pos) > 0){
pos.interaction <- paste(mental.sheet[i, Box], '->',
names(mental.sheet[i, pos + 1, with = F]),
sep = '')
model <- append(model, pos.interaction)
}
neg <- which(mental.sheet[i, .SD, .SDcol = box.names] < 0)
if(length(neg) > 0){
neg.interaction <- paste(mental.sheet[i, Box], '-*',
names(mental.sheet[i, neg + 1, with = F]),
sep = '')
model <- append(model, neg.interaction)
}
}
return(model)
}
#function to create qnm models for qpress from signed digraphs
make.qnm<-function(modmat){
q<-MM2Qpress(modmat)
qnm<-dput(q)
qnm.mod<-parse.digraph(qnm)
qnm.model<-enforce.limitation(qnm.mod)
}
#--------------------------------------------------------------------
#####PERTURBATIONS#####
#to loop over multiple models
#create folder with files in directory, call list of files
set.seed(1984)
#get(load RData files if required
#create list
# mod<-list(WSS62qnm10, WSS62qnm20, WSS62qnm30, WSS62qnm40, WSS62qnm50, WSS62qnm100)
mod<-list(WSS62qnm50)
#get MM matrix into qpress format, create model
q.models<-vector('list', length(mod))
for(i in seq_along(mod)){
q.models[[i]]<-make.qnm(mod[[i]])
}
# #write dia files
# for(i in seq_along(q.models)){
# #can add text below to further specify name
# write.dia(q.models[[i]], paste0("WSS", i,".dia"))
# }
#loop through models to create results csv file and pdf barplots
for(j in seq_along(q.models)){
A <- adjacency.matrix(q.models[[j]])
#A
## Function to generate the community matrix
s <- community.sampler(q.models[[j]])
## Function to check the validation condition
#press <- press.validate(GB2el,
# perturb=c("Forage_Fish"=-1),
# monitor=c("Commercial_Pelagic_Fishery"=1))
#Press/Perturbation with decrease forage fish and no monitoring)
press <- press.validate(q.models[[j]],
perturb=c("Phytoplankton"=1), monitor=F)
## Function to define the perturbation scenario
impact <- press.impact(q.models[[j]],perturb=c("Phytoplankton"=1))
## Use 1000 simulations
n.sims <- 1000
results <- 0
i <- 0
while(i < n.sims) {
## Randomly choose edges to retain
z <- s$select(runif(1))
## Sample community matrix
W <- s$community()
## Check press condition and stability
if(!(press(W) && stable.community(W))) next
## Monitor impact post press
imp <- impact(W)
imp[abs(imp)<1e-07]=0 #added to remove small perturbations
results <- results + outer(sign(imp),-1:1,'==')
i <- i+1
rownames(results) <- levels(q.models[[j]]$From)
colnames(results) <- c('-','0','+')
# save(results, paste0("WSS62","results", j,"Phytoplankton_up",".RData"))
write.csv(results, paste0("WSS62","results", j,"Phytoplankton_up",".csv"))
library(RColorBrewer)
pal <- brewer.pal(n=5,"Greys")[5:1]
#opar <- par(mar=c(5,10,1,1)+0.1)
prop <- results/rowSums(results)
r <- colSums(t(prop)*(-1:1))
filename= paste("WSS62", "plots", j,"Phytoplankton_up", ".pdf")
pdf(filename)
#adjust barplot if you want it ordered by proportion
par(mar=c(5,10,1,1)+0.1)
barplot(t(prop),
horiz=T,cex.names=0.7,cex.axis=0.8,las=2,border=T, col=pal,xlab="Proportion")
dev.off()
}
}
#code for ordered barplot
# pal <- brewer.pal(n=5,"Greys")[5:1]
# opar <- par(mar=c(5,10,1,1)+0.1)
# prop <- results/rowSums(results)
# r <- colSums(t(prop)*(-1:1))
# filename= paste("wss", "plot", j,"phytoplankton up", ".pdf")
# pdf(filename)
# barplot(t(prop[order(r),]),
# horiz=T,cex.names=0.6,cex.axis=0.8,las=2,border=T,col=pal,xlab="Proportion")
# par(opar)
gather results
#gather positive simulation results
WSS28.results_pos<-WSS28.results %>%
filter(Outcome=="Positive") %>%
unite("model_ID", Group:Direction, remove=FALSE) %>%
select(!2:5)%>%
pivot_longer(!model_ID, names_to = "species", values_to = "pos")
#gather neutral simulation results, calculate half
WSS28.results_neu<-WSS28.results %>%
filter(Outcome=="Neutral")%>%
unite("model_ID", Group:Direction, remove=FALSE) %>%
select(!2:5)%>%
pivot_longer(!model_ID, names_to = "species", values_to = "neu") %>%
mutate(half_neu=neu/2)
#join positive and neutral simulations
WSS28.results_posneu<-left_join(WSS28.results_pos, WSS28.results_neu) %>%
mutate(result=pos+half_neu) %>%
# mutate(result=(positive+500)/2) %>%
select(model_ID, species,result)%>%
pivot_wider(names_from=species, values_from=result, names_vary='slowest')
#create columns that match with QNM results
WSS28.cols<-WSS28.results %>%
filter(Outcome=="Positive") %>%
select(1:2)
WSS28.results.2<-bind_cols(WSS28.cols, WSS28.results_posneu)
#align results, correct for differences in spelling
WSS28.results_ecos<-WSS28.results.2 %>%
mutate(
Direction=case_match(Direction,"Plus" ~ "plus", "Minus" ~ "minus"),
Group=case_match(Group, "Seals" ~ "Seal", .default=Group)
) %>%
add_column (Model="ecosense", model_links=0) %>%
select(!model_ID)%>%
unite("model_ID", Group:Direction, remove=FALSE)
# WSS28.results_ecos<-WSS28.results %>%
# mutate(
# Direction=case_match(Direction,"Plus" ~ "plus", "Minus" ~ "minus"),
# Group=case_match(Group, "Seals" ~ "Seal", .default=Group)
# ) %>%
# select(!Outside)%>%
# filter(Outcome=="Positive") %>%
# add_column (Model="ecosense", model_links=0) %>%
# unite("model_ID", Group:Direction, remove=FALSE)
#Bind Ecosense and QNM results
WSS28_results_ecosQNM<-bind_rows(WSS28.results_ecos, WSS28_QNM_results_compiled_tr)
saveRDS(WSS28_results_ecosQNM, "~/JCT work/WGNARS/RpathtoQNM/RpathQNM/data/WSS28_results_ecosQNM_wfishery.RDS")
Graphing ecosense and qnm (THIS IS OUTDATED)
library(reshape2)
library(ggplot2)
library(gridExtra)
require(tidyverse)
WSS28_results_ecosQNM <- readRDS("~/JCT work/WGNARS/RpathtoQNM/RpathQNM/data/WSS28_results_ecosQNM.RDS")
ecoQNMresults<-WSS28_results_ecosQNM %>%
unite("scenario", Group:Direction, remove=FALSE ) %>%
unite("Model_name", Model:model_links, remove=FALSE)
ecoQNMresults_long<-ecoQNMresults %>%
pivot_longer (cols= 8:35, names_to="node" ) %>%
mutate(sign=value) %>%
mutate(sign=case_when(sign >= 600 ~ "positive",
sign <600 & sign >=400~ "zero",
sign <400 ~ "negative")) %>%
# mutate(fill=value) %>%
mutate(scaled_value= (abs(value-500)))
#
# ecoQNMresults_long$fill[ecoQNMresults_long$fill >800] <- "strong"
# ecoQNMresults_long$fill[ecoQNMresults_long$fill >600 & ecoQNMresults_long$fill <=800] <- "weak"
# ecoQNMresults_long$fill[ecoQNMresults_long$fill >400 & ecoQNMresults_long$fill <=600] <- "weak"
# ecoQNMresults_long$fill[ecoQNMresults_long$fill >200 & ecoQNMresults_long$fill <=400] <- "weak"
# ecoQNMresults_long$fill[ecoQNMresults_long$fill <=200] <- "strong"
########################
# Plot up responses as bubble plots
dat_fishplus<- ecoQNMresults_long %>%
filter(scenario=="Fishery_plus")
p_fishplus<- ggplot(dat_fishplus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black", "white"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black", "white")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
dat_fishminus<- ecoQNMresults_long %>%
filter(scenario=="Fishery_minus")
p_fishminus<- ggplot(dat_fishminus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
dat_sealminus<- ecoQNMresults_long %>%
filter(scenario=="Seal_minus")
p_sealminus<- ggplot(dat_sealminus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
dat_sealplus<- ecoQNMresults_long %>%
filter(scenario=="Seal_plus")
p_sealplus<- ggplot(dat_sealplus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
dat_smallpelagicsminus<- ecoQNMresults_long %>%
filter(scenario=="Small pelagics_minus")
p_smallpelagicsminus<- ggplot(dat_smallpelagicsminus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
dat_smallpelagicsplus<- ecoQNMresults_long %>%
filter(scenario=="Small pelagics_plus")
p_smallpelagicsplus<- ggplot(dat_smallpelagicsplus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
dat_phytoplanktonminus<- ecoQNMresults_long %>%
filter(scenario=="Phytoplankton_minus")
p_phytoplanktonminus<- ggplot(dat_phytoplanktonminus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
dat_phytoplanktonplus<- ecoQNMresults_long %>%
filter(scenario=="Phytoplankton_plus")
p_phytoplanktonplus<- ggplot(dat_phytoplanktonplus)+
geom_point(aes(y=node, x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0)
)
-------------------------------------------
dat_fishminus<- ecoQNMresults_long %>%
filter(scenario=="Fishery_minus")
p_fishminus<- ggplot(dat_fishminus)+
geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25, 24)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
labs(title="fishery -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_sealminus<- ecoQNMresults_long %>%
filter(scenario=="Seal_minus")
p_sealminus<- ggplot(dat_sealminus)+
geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","black"))+
scale_shape_manual(values=c(25,4)) +
scale_fill_manual(values=c("paleturquoise4", "black")) +
labs(title="seal -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_sealplus<- ecoQNMresults_long %>%
filter(scenario=="Seal_plus")
p_sealplus<- ggplot(dat_sealplus)+
geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
labs(title="seal +")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_smallpelagicsminus<- ecoQNMresults_long %>%
filter(scenario=="Small pelagics_minus")
p_smallpelagicsminus<- ggplot(dat_smallpelagicsminus)+
geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
labs(title="small pelagic -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_smallpelagicsplus<- ecoQNMresults_long %>%
filter(scenario=="Small pelagics_plus")
p_smallpelagicsplus<- ggplot(dat_smallpelagicsplus)+
geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
labs(title="small pelagic +")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_phytoplanktonminus<- ecoQNMresults_long %>%
filter(scenario=="Phytoplankton_minus")
p_phytoplanktonminus<- ggplot(dat_phytoplanktonminus)+
geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
labs(title="phytoplankton -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_phytoplanktonplus<- ecoQNMresults_long %>%
filter(scenario=="Phytoplankton_plus")
p_phytoplanktonplus<- ggplot(dat_phytoplanktonplus)+
geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value,
pch=sign, color=sign, fill=sign)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+
scale_shape_manual(values=c(25,24,4)) +
scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) +
labs(title="phytoplankton +")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
Reorganizes the combined tables with Rpath and QNM results
library(data.table); library(here)
#QNM results
files <- as.data.table(dir(here::here('WSS28model', 'data neutrals')))
#only use .RData
files <- files[V1 %like% '.RData']
all.results <- c()
for(isim in 1:nrow(files)){
load(here::here('WSS28model', 'data neutrals', files[isim, ]))
#Create output data.tablereg
model <- gsub(".*neu(.*)\\..*", "\\1", files[isim, ])
model_ID <- paste('QNM_', substr(model, 1, 1), 0, sep = '')
scenario <- substr(model, 2, 100)
groups <- dimnames(results)[[1]]
neg <- results[, 1]
neu <- results[, 2]
pos <- results[, 3]
out <- data.table(Model_ID = model_ID,
Scenario = scenario,
Group = groups,
Neg = neg,
Neu = neu,
Pos = pos)
#Assign Negative, Neutral, Positive, or Mixed result
out[Neg > 600, Symbol := 'Neg']
out[Neu > 600, Symbol := 'Neu']
out[Pos > 600, Symbol := 'Pos']
out[is.na(Symbol), Symbol := 'Mix']
#Assign Strength of dominant direction
out[Symbol == 'Neg', Size := Neg + 0.5 * Neu]
out[Symbol == 'Pos', Size := Pos + 0.5 * Neu]
out[Symbol == 'Neu', Size := 701]
out[Symbol == 'Mix', Size := 701]
out[Size > 599 & Size <= 700, Strength := 'Weak']
out[Size > 700 & Size <= 900, Strength := 'Moderate']
out[Size > 900, Strength := 'Strong']
#Drop extra columns
#out[, c('Neg', 'Neu', 'Pos') := NULL]
#join
all.results <- rbindlist(list(all.results, out))
}
#Change QNM_60 to QNM_0
all.results[Model_ID == 'QNM_60', Model_ID := 'QNM_0']
#Ecosense results
load(here::here('data', 'WSS28.results.rda'))
groups <- colnames(WSS28.results)[which(!colnames(WSS28.results) %in%
c('Group', 'Direction', 'Outcome', 'Outside'))]
#Clean up names to match QNM
WSS28.results[Direction == 'Plus', Direction := 'up']
WSS28.results[Direction == 'Minus', Direction := 'down']
WSS28.results[Group == 'Small pelagics', Group := 'Spelagics']
for(isp in 1:length(groups)){
sp.result <- data.table::dcast(WSS28.results, Group + Direction ~ Outcome,
value.var = groups[isp])
sp.result[, Model_ID := 'Rpath']
sp.result[, Scenario := paste(Group, Direction, sep = '_')]
sp.result[, Group := NULL]
sp.result[, Group := groups[isp]]
#Assign Negative, Neutral, Positive, or Mixed result
sp.result[Negative > 600, Symbol := 'Neg']
sp.result[Neutral > 600, Symbol := 'Neu']
sp.result[Positive > 600, Symbol := 'Pos']
sp.result[is.na(Symbol), Symbol := 'Mix']
#Assign Strength of dominant direction
sp.result[Symbol == 'Neg', Size := Negative + 0.5 * Neutral]
sp.result[Symbol == 'Pos', Size := Positive + 0.5 * Neutral]
sp.result[Symbol == 'Neu', Size := 701]
sp.result[Symbol == 'Mix', Size := 701]
sp.result[Size > 599 & Size <= 700, Strength := 'Weak']
sp.result[Size > 700 & Size <= 900, Strength := 'Moderate']
sp.result[Size > 900, Strength := 'Strong']
#Drop extra columns
setnames(sp.result, c('Negative', 'Neutral', 'Positive'), c('Neg', 'Neu', 'Pos'))
sp.result[, 'Direction' := NULL]
#join
all.results <- rbindlist(list(all.results, sp.result), use.names = T)
}
data.table::setkey(all.results, 'Model_ID', 'Scenario', 'Group')
usethis::use_data(all.results, overwrite = T)
Graphing combined plots
library(reshape2)
library(ggplot2)
library(gridExtra)
require(tidyverse)
require(grid)
library(here)
require(stringr)
load("~/JCT work/WGNARS/RpathtoQNM/RpathQNM/data/all.results.rda")
group_order<-c("Detritus", "Discards", "Phytoplankton", "Microflora", "Microzoop", "Mesozoop", "Macrozoop", "Gelatinous zoop", "Meiofauna", "Worms", "Other molluscs", "Bivalves", "Shrimps", "Small crab other arthropoda", "Megabenthos", "Squids", "Other pelagic mesopelagic", "Small pelagics", "Skates", "L benthivores", "D benthivores", "D piscivores", "Large pelagics", "Sharks", "Sea birds", "Seals", "Toothed cetaceans", "Whales", "Fishery")
model_order<-c("Rpath","QNM_0",
"QNM_10",
"QNM_20",
"QNM_30",
"QNM_40",
"QNM_50")
all.results_noneut<-as_tibble(all.results)
#change the naming of Model_ID, this has been fixed in table code so eleminate if need be.
all.results_noneut<-all.results_noneut%>%
mutate_if(is.character, str_replace_all, pattern="QNM_1", replacement="QNM_10") %>%
mutate_if(is.character, str_replace_all, pattern="QNM_2", replacement="QNM_20") %>%
mutate_if(is.character, str_replace_all, pattern="QNM_3", replacement="QNM_30") %>%
mutate_if(is.character, str_replace_all, pattern="QNM_4", replacement="QNM_40") %>%
mutate_if(is.character, str_replace_all, pattern="QNM_5", replacement="QNM_50") %>%
mutate_if(is.character, str_replace_all, pattern="QNM_6", replacement="QNM_0")
# Plot up responses as bubble plots
dat_fishplus<- all.results_noneut %>%
filter(Scenario=="Fishery_up")
dat_fishplus<-dat_fishplus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_fishplus<- ggplot(dat_fishplus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Fishery +")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_fishminus<- all.results_noneut %>%
filter(Scenario=="Fishery_down")
dat_fishminus<-dat_fishminus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_fishminus<- ggplot(dat_fishminus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Fishery -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_sealplus<- all.results_noneut %>%
filter(Scenario=="Seals_up")
dat_sealplus<-dat_sealplus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_sealplus<- ggplot(dat_sealplus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Seal +")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_sealminus<- all.results_noneut %>%
filter(Scenario=="Seals_down")
dat_sealminus<-dat_sealminus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_sealminus<- ggplot(dat_sealminus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Seals -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_smallpelagicsplus<- all.results_noneut %>%
filter(Scenario=="Spelagics_up")
dat_smallpelagicsplus<-dat_smallpelagicsplus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_smallpelagicsplus<- ggplot(dat_smallpelagicsplus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Small pelagics +")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_smallpelagicsminus<- all.results_noneut %>%
filter(Scenario=="Spelagics_down")
dat_smallpelagicsminus<-dat_smallpelagicsminus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_smallpelagicsminus<- ggplot(dat_smallpelagicsminus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Small pelagics -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_phytoplanktonplus<- all.results_noneut %>%
filter(Scenario=="Phytoplankton_up")
dat_phytoplanktonplus<-dat_phytoplanktonplus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_phytoplanktonplus<- ggplot(dat_phytoplanktonplus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Phytoplankton +")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
dat_phytoplanktonminus<- all.results_noneut %>%
filter(Scenario=="Phytoplankton_down")
dat_phytoplanktonminus<-dat_phytoplanktonminus %>%
arrange(match(Group, group_order))%>%
arrange(match(Model_ID, model_order))
p_phytoplanktonminus<- ggplot(dat_phytoplanktonminus)+
geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size,
pch=Symbol, color=Symbol, fill=Symbol)) +
# pch=sign, color=sign, fill=fill)) +
scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+
scale_shape_manual(values=c(4, 25, 20,24)) +
scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) +
labs(title="Phytoplankton -")+
theme_bw() +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.text.x=element_text (angle=-90, hjust=0),
legend.position="none"
)
grid.newpage()
jpeg(filename="RpathQNM_combined_plots_revisedRpath_2025.jpeg", width=20, height=15, units="in", quality=100,res=600)
grid.arrange(p_phytoplanktonplus, p_smallpelagicsplus, p_sealplus, p_fishplus, p_phytoplanktonminus, p_smallpelagicsminus, p_sealminus, p_fishminus, ncol = 4)
calculate outdegrees
#for perturbed nodes in all QNM models
library(tidyverse); library(here)
Model_ID<-c("Rpath", "QNM_0", "QNM_10", "QNM_20", "QNM_30", "QNM_40", "QNM_50")
outd_qnm_10<-(WSS28qnm10[,-1])
rownames(outd_qnm_10)<-WSS28qnm10$FG
outd_qnm_10[outd_qnm_10>0]<-1
outd_qnm_10[outd_qnm_10<0]<-0
outd_qnm_20<-(t(SS28qnm20))
outd_qnm_20[outd_qnm_20>0]<-1
outd_qnm_20[outd_qnm_20<0]<-0
outd_qnm_30<-WSS28qnm30
outd_qnm_30[outd_qnm_30>0]<-1
outd_qnm_30[outd_qnm_30<0]<-0
outd_qnm_40<-WSS28qnm40
outd_qnm_40[outd_qnm_40>0]<-1
outd_qnm_40[outd_qnm_40<0]<-0
outd_qnm_50<-WSS28qnm50
outd_qnm_50[outd_qnm_50>0]<-1
outd_qnm_50[outd_qnm_50<0]<-0
outd_qnm_0<-WSS28qnm100
outd_qnm_0[outd_qnm_0>0]<-1
outd_qnm_0[outd_qnm_0<0]<-0
outd_QNM_10p<-outd_qnm_10%>% filter ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum))
outd_QNM_20p<-outd_qnm_20%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum))
outd_QNM_30p<-outd_qnm_30%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum))
outd_QNM_40p<-outd_qnm_40%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum))
outd_QNM_50p<-outd_qnm_50%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum))
outd_QNM_0p<-outd_qnm_0%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum))
ind_rpathQNM<-rbind(outd_QNM_0p, outd_QNM_0p, outd_QNM_10p, outd_QNM_20p, outd_QNM_30p, outd_QNM_40p, outd_QNM_50p)
ind_rpathQNM<-cbind(ind_rpathQNM, Model_ID)
degrees<-rep(c("indegree", "outdegree", "in_out_degree"), each=7)
alldeg_rpathQNM<-rbind(ind_rpathQNM, outd_rpathQNM, deg_rpathQNM)
alldeg_rpathQNM<-cbind(alldeg_rpathQNM, degrees)
save(alldeg_rpathQNM, file="data/alldegrees_rpathQNM.RData")
NOAA-EDAB/RpathQNM documentation built on June 13, 2025, 3:02 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Setup
require(devtools) require(Rtools) require(data.table) require (Rpath) require(QPress) require(tidyverse) # install_github("NOAA-EDAB/Rpath") # install_github("SWotherspoon/QPress")
Rpath for WSS28
#Western Scotian Shelf Model (No stanza groups) #SML #User parameters #------------------------------------------------------------------------------- #Required packages library(data.table); library(Rpath) setwd("~/JCT work/WGNARS/RpathtoQNM/RpathQNM/WSS28model") #------------------------------------------------------------------------------- #User created functions #------------------------------------------------------------------------------- #Western Scotian Shelf Groups - collapsed list (no multistanza groups) groups <- c('Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks', 'Large pelagics', 'D piscivores', 'L benthivores', 'Skates', 'Small pelagics', 'Other pelagic mesopelagic', 'Squids', 'Megabenthos', 'Small crab other arthropoda', 'Shrimps', 'Bivalves', 'Other molluscs', 'Worms', 'Meiofauna','Gelatinous zoop', 'Macrozoop', 'Mesozoop', 'Microzoop', 'Microflora', 'Phytoplankton', 'Discards', 'Detritus', 'Fishery') #Identify type (0 = consumer, 1 = producer, 2 = detritus, 3 = fleet) types <- c(rep(0, 24), 1, 2, 2, 3) #Create parameter list object WSS28.params <- create.rpath.params(groups, types) #Fill in model parameters #Biomass, Production, consumption biomass <- c(0.407075, 0.049754, 0.04383, 0.00617, 0.0260249,0.023573, 4.082286, 4.240253, 0.137757, 6.294171, 0.878285, 0.173893, 0.458598, 1.818481, 1.302563, 64.01978, 2, 7.345485, 42.75382, 0.5204, 41.31, 23.79555, 5.800508, 3.678861, 33.664, 0.063293, 1,NA) pb <- c(0.071,0.18,0.147007,0.25,0.18,0.4,0.5341109,0.3798308,0.2320817,0.6959251,0.7432774,4, 0.8221899,2.759435,3,0.69,0.75,1.25,1.32583,15.51,3.04,29.2, 82.8, 104.938, 70.639, rep(NA, 3)) qb <- c(4.94,14.5,7.338632,87.6,4.78,4.24,2.868397,3.254083,2.447539,3.475426,3.532147,11.33333, 5.481265,18.39623 ,rep(NA, 4), 8.83887,62.05, 19.5, 73, rep(NA, 6)) WSS28.params$model[, Biomass := biomass] WSS28.params$model[, PB := pb] WSS28.params$model[, QB := qb] #Production to Consumption WSS28.params$model[Group %in% c('megabenthos', 'Shrimps', 'Bivalves', 'Other molluscs', 'Worms' ), ProdCons := 0.15] #WSS28.params$model[Group == 'Mesozoop', ProdCons := 0.4] WSS28.params$model[Group %in% c('Microzoop', 'Microflora'), ProdCons := 0.5] #Biomass accumulation and unassimilated production ba <- c(0.0032566,0,0.003702101,0,-0.001951868,0,-0.05736,0.044246, rep(0, 4), 0.0301295, 0,0, 0.03234285,rep(0,11), NA) # ba<-c(rep(0,28), NA) WSS28.params$model[, BioAcc := ba] WSS28.params$model[, Unassim := c(rep(0.2, 20), rep(0.4, 3), 0.2, rep(0, 3), NA)] #Detrital fate WSS28.params$model[, Detritus := c(rep(1, 25), rep(0, 3))] WSS28.params$model[, Discards := c(rep(0, 27), 1)] #Landings/Discards land <- c(rep(0, 4),0.00209077, 0.005067847, 0.3251713, 0.182276, 0.000799, 0.999061, 0.019708, 0.001637, 0.220424,0, 0.000162, 0.117151, 0.00236,0, 0.03026,rep(0,8),NA) disc <- c(0,0.002995,0,2.4E-05, 0.000575, 0.000469, 0.013477, 0.002266, 0.013094, 0.000487, 0.00056, 1.9E-05, 0.008776, 0.000633, 8E-06 ,0,4.1E-05,0,0.019869, rep(0, 8), NA) WSS28.params$model[, Fishery := land] WSS28.params$model[, Fishery.disc := disc] #Diet composition DC.groups <- groups[which(!groups %in% c('Discards', 'Detritus', 'Fishery'))] # DC.data <- as.data.table(matrix(c(0, 0, 0, 0 ,0, 0 ,0, 0, 0, 0, 0,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0, 0,0,0,0,0,0, # 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0, # 0, 0 ,0, 0 ,0.006511993, 0 ,0, 0 ,0, 0 ,0, 0, 0, 0 ,0, 0 ,0, 0, 0, 0 ,0, 0, 0, 0, 0, # 0, 0, 0, 0.002663992, 0, 0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0 ,0 ,0, 0 ,0, 0 ,0, 0, 0, # 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0 ,0, 0, # 0, 0, 0, 0 ,0, 0.002219998 ,0, 0, 0, 0, 0, 0 ,0, 0, 0, 0 ,0, 0 ,0, 0 ,0, 0, 0, 0, 0, # 0.004802, 0.2437588, 0.1276689, 0.0007209978, 0.1848318, 0.06917194, 0.05424953, 0.005541045 ,0.003080484, 0 ,0, 0.003657, 0, 0, 0, 0 ,0 ,0, 0, 0,0, 0, 0, 0 ,0, # 0.012545, 0.09793591, 0.2855617, 0.005076986, 0.3280216, 0.1511799, 0.04674805, 0.004132936, 0.05460171, 0, 0, 0.007678, 0.001092015, 0, 0, 0 ,0, 0 ,0 ,0, 0, 0 ,0, 0, 0, # 0, 0 ,0.006431994, 0, 0.0008079991, 0.002838997, 0.0004965155, 0 ,0, 0 ,0, 0.000187, 0, 0 ,0, 0, 0, 0 ,0, 0, 0, 0, 0, 0, 0, # 0.173918, 0.3776346 ,0.5207465, 0.3529089, 0.4223085, 0.6149244, 0.1726697, 0.00143368, 0.02304704, 0, 0, 0.015327, 0, 0, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0, 0, 0, # 0.0194, 0.007396993, 0.05680194, 0.07901775, 0.04387695, 0.1276559, 0.02289399, 0.01478033, 0.03303853, 0, 0, 0.004734, 0 ,0, 0, 0 ,0 ,0 ,0, 0, 0 ,0, 0 ,0, 0, # 0.011769, 0.2679567, 0.002056998, 0.03595289, 0.01354798, 0.01555398, 0.02838177, 0.001817405 ,0.006008342, 0, 0, 0.037929, 0.002280137, 0, 0, 0, 0, 0 ,0, 0, 0, 0 ,0, 0, 0, # 0, 0 ,2.099998E-05, 0, 9.29999E-05, 7.299994E-05, 0.003771638, 0.001804974, 0.00659312, 0, 0, 0, 0.007052792, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # 0, 0, 0, 0, 0, 0, 0.03455577, 0.07374978, 0.1577646, 0.02415398, 0.060237, 0.014737, 0.09158555, 0.009690308, 0.004533, 0, 0 ,0.03769, 0.02, 0.0001840747, 0, 0, 0, 0, 0, # 0, 0, 0.0007109993 ,0, 0, 0.003309997, 0.07398107, 0.08974432, 0.07314944, 0.0622156, 0.03105934, 0.010014, 0.001267644, 0 ,0 ,0 ,0 ,0.020851, 0, 0, 0, 0, 0, 0, 0, # 0, 0 ,0 ,0 ,0 ,0 ,0.0008845589, 0.00388914, 0, 0, 0, 0, 0.004711973, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0 ,0 ,0 ,0, # 0, 0, 0, 0, 0, 0 ,0.0006980912, 0.02864456, 0.0008348122, 0 ,0, 0 ,0.3265275 ,0.01410418 ,0 ,0 ,0, 0.332281, 0, 0.001184787, 0, 0 ,0 ,0 ,0, # 0, 0, 0, 0, 0, 0, 0.001044587, 0.005987847, 0.002030869 ,0, 0 ,0 ,0.0419751, 0.01410418, 0 ,0 ,0 ,0.027689 ,0, 0.0006171872, 0 ,0, 0, 0, 0, # 0, 0, 0 ,0 ,0, 0, 0.002290674, 0.09854721, 0.1947723, 0.001041604, 0.0008129025, 0.013273, 0.1380028, 0.02820811, 0.015272, 0 ,0 ,0, 0.070001, 0.0006171872, 0, 0, 0, 0, 0, # 0, 0, 0, 0, 0, 0.01216899, 0.004603415, 0.09481322, 0.001261869, 0.003441011, 0, 0, 0.281883, 0.08955499, 0 ,0 ,0,0.304589, 0.070001, 0.03959263, 0, 0, 0, 0, 0, # 0, 0, 0, 0, 0, 5.899995E-05, 0.08581191, 0.008719152, 0.001671519, 0.0001719213, 0.001211964, 0.001571, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # 0.726929, 0, 0, 0.315373, 0, 0, 0.4453562, 0.5160776, 0.3443648, 0.8572081, 0.3833565, 0.890893, 0, 0 ,0.123605, 0 ,0 ,0 ,0 ,0.005774691, 0.161448, 0.055941, 0, 0 ,0, # 0.050637, 0 ,0, 0, 0, 0, 0.01242108, 0.0372947, 0.09778056, 0.05176783, 0.5141716, 0, 0, 0, 0.244326, 0 ,0 ,0 ,0 ,0.005284418 ,0.696038 ,0.47728, 0.058572 ,0, 0, # 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.00386473, 0.095929, 0.11221, 0.161911, 0.099731, 0, # 0, 0, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0 ,0 ,0 ,0.006854988, 0.008578, 0.04999995, 0.095, 0, 0, 0.0102941, 0.005032, 0.015172, 0.076645, 0.202492, 0, # 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.002426699, 0, 0, 0, 0.069225, 0.4222216, 0.882656, 0 ,0, 0.08967555, 0.041553, 0.127274, 0.702872, 0.209335, 0, # 0 ,0.005316995, 0, 0.0005529983, 0, 0, 0, 0, 0, 0, 0, 0, 0.002525342, 0.001344179, 0, 0, 0, 0, 0, 0 ,0 ,0 ,0, 0, 0, # 0 ,0 ,0, 0 ,0, 0, 0, 0, 0, 0 ,0, 0, 0.06347585, 0.8361391, 0.534461, 0.5277785, 0.022344, 0.2769, 0.839998, 0.8429105, 0, 0.212123, 0, 0.488442, 1, # 0, 0, 0, 0.2077324, 0, 0.000843999, 0.009141359, 0.01302205, 0, 0, 0.006723996, 0, 0.03762027, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), # 29,25, byrow = T)) #---------------------------- # code to copy paste data into Rpath, can do this or copy paste to CSV and import. # read.excel <- function(header=TRUE,...) { # read.table("clipboard",sep="\t",header=F,...) # } # # dcdat<-read.excel() # DC.data<-as.data.table(dcdat) # save(DC.data, file='DC.data.noscallop.RData') load('DC.data.noscallop.RData') #------------------------------------------- #Fix column names # setnames(WSS28.params$diet, 1, 'Group') # setnames(WSS28.params$diet, paste0('V', 1:25), DC.groups) # setnames(WSS28.params$diet, c(paste('V', 1:27, sep = '')), c('Group', DC.groups) #Merge DC matrix with diet parameter object WSS28.params$diet <- cbind(WSS28.params$diet[, Group], DC.data) #Fix column names setnames(WSS28.params$diet, 1, 'Group') setnames(WSS28.params$diet, paste0('V', 1:25), DC.groups) #Ecopath WSS28 <- rpath(WSS28.params, 'Western Scotian Shelf') save(WSS28.params, file = "WSS28.RData")
Make CModels WSS28
#QPress------------------------------------------------------------------------- setwd("~/JCT work/WGNARS/RpathQNM/RpathQNM/WSS28model") #Prey matrix qpress.prey <- WSS28.params$diet pred.cols <- names(qpress.prey)[which(names(qpress.prey) != 'Group')] qpress.prey <- qpress.prey[Group != 'Import'] #Ensure that the matrices are square detritus.DC <- data.table(Group = qpress.prey[, Group], Detritus = 0, Discards = 0) qpress.prey <- merge(qpress.prey, detritus.DC, by = 'Group') #line columns and rows up groupnames <- as.character(qpress.prey$Group) groupnames<-c('Group', groupnames) setcolorder(qpress.prey, groupnames) #Predator M2/F matrix x <- copy(WSS28) ngroup <- x$NUM_LIVING + x$NUM_DEAD out <- data.frame(Group = x$Group[1:ngroup], type = x$type [1:ngroup], PB = x$PB [1:ngroup]) #Calculate F mortality # totcatch <- x$Catch + x$Discards totcatch <- x$Landings + x$Discards Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))]) # setnames(Fmort, 'Fleet1', 'Fishery') # Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))]) # setnames(Fmort, paste('V', 1:x$NUM_GEARS, sep = ''), # paste('F.', x$Group[(ngroup +1):x$NUM_GROUPS], sep = '')) out <- cbind(out, Fmort[1:ngroup, ]) setnames(out, 'Fmort[1:ngroup, ]', 'Fishery') # pred.cols <- names(qpress.prey)[which(names(qpress.prey) != 'Group')] # # qpress.prey <- qpress.prey[Group != 'Import'] # # #Ensure that the matrices are square # detritus.DC <- data.table(Group = qpress.prey[, Group], Detritus = 0, Discards = 0) # qpress.prey <- merge(qpress.prey, detritus.DC, by = 'Group') # # #line columns and rows up # groupnames <- as.character(qpress.prey$Group) # groupnames<-c('Group', groupnames) # setcolorder(qpress.prey, groupnames) # # #Predator M2/F matrix # x <- copy(WSS28) # ngroup <- x$NUM_LIVING + x$NUM_DEADs # out <- data.frame(Group = x$Group[1:ngroup], # type = x$type [1:ngroup], # PB = x$PB [1:ngroup]) # #Calculate F mortality # totcatch <- x$Catch + x$Discards # Fmort <- as.data.frame(totcatch / x$BB[row(as.matrix(totcatch))]) # setnames(Fmort, paste('V', 1:x$NUM_GEARS, sep = ''), # paste('F.', x$Group[(ngroup +1):x$NUM_GROUPS], sep = '')) # out <- cbind(out, Fmort[1:ngroup, ]) # setnames(out, 'Fmort[1:ngroup, ]', 'Fishery') #Add Fishery to prey fishery <- as.data.table(out) fishery <- fishery[, list(Group, Fishery)] qpress.prey <- merge(qpress.prey, fishery, by="Group") qpress.prey<-qpress.prey[, -1]#remove Group column fishery.0 <- as.data.table(matrix(rep(0, nrow(fishery) + 1), 1, nrow(fishery) + 1)) # setnames(fishery.0, paste('V', 1:x$NUM_Groups, sep=''), groupnames) qpress.prey <- rbindlist(list(qpress.prey, fishery.0), use.names=FALSE) # fishery <- as.data.table(out) # fishery <- fishery[, list(Group, Fishery)] # qpress.prey <- merge(qpress.prey, fishery, by = 'Group') # qpress.prey[, Group := NULL] # fishery.0 <- as.data.table(matrix(rep(0, nrow(fishery) + 1), 1, nrow(fishery) + 1)) # qpress.prey <- rbindlist(list(qpress.prey, fishery.0)) # #Predator M2/F matrix # x <- copy(WSS28.b) # ngroup <- x$NUM_LIVING + x$NUM_DEAD # out <- data.frame(Group = x$Group[1:ngroup], # type = x$type [1:ngroup], # PB = x$PB [1:ngroup]) #Calculate M0 M0 <- c(x$PB[1:x$NUM_LIVING] * (1 - x$EE[1:x$NUM_LIVING]), x$EE[(x$NUM_LIVING + 1):ngroup]) out <- cbind(out, M0) #Calculate M2 bio <- x$Biomass[1:x$NUM_LIVING] BQB <- bio * x$QB[1:x$NUM_LIVING] diet <- as.data.frame(x$DC) nodetrdiet <- diet[1:x$NUM_LIVING, ] detrdiet <- diet[(x$NUM_LIVING +1):ngroup, ] newcons <- nodetrdiet * BQB[col(as.matrix(nodetrdiet))] predM <- newcons / bio[row(as.matrix(newcons))] detcons <- detrdiet * BQB[col(as.matrix(detrdiet))] predM <- rbind(predM, detcons) # setnames(predM, paste('V', 1:x$NUM_LIVING, sep = ''), # paste('M2.', x$Group[1:x$NUM_LIVING], sep = '')) qpress.pred <- as.data.table(cbind(out, predM)) M2.cols <- names(qpress.pred)[which(!names(qpress.pred) %in% c('Group', 'type', 'PB', 'M0'))] qpress.pred[, mort.sum := rowSums(.SD), .SDcols = M2.cols] qpress.pred <- qpress.pred[, .SD / mort.sum, .SDcols = M2.cols, by = Group] qpress.pred[is.na(qpress.pred)] <- 0 #Code for thresholds # for(ipred in seq_along(M2.cols)){ # setnames(qpress.pred, M2.cols[ipred], 'V1') # qpress.pred[V1 < .1, V1 := 0] # setnames(qpress.pred, 'V1', M2.cols[ipred]) # } #make qpress.pred negative qpress.pred <- qpress.pred[, .SD * -1, .SDcol = M2.cols, by = Group] #Ensure that the matrices are square detritus.m2 <- data.table(qpress.pred[, 'Group'], M2.Detritus = 0, M2.Discards = 0) qpress.pred <- merge(qpress.pred, detritus.m2, by = 'Group') qpress.pred <- qpress.pred[, .SD, .SDcols = c(M2.cols, 'M2.Detritus', 'M2.Discards')] prednames <- as.character(c('Fishery', detritus.DC$Group)) colnames(qpress.pred) <- prednames new.names<-colnames(qpress.prey) setcolorder(qpress.pred, new.names) qpress.pred <- t(qpress.pred) qpress.pred <- cbind(qpress.pred, rep(0, nrow(qpress.pred))) colnames(qpress.pred) <- c(new.names) #Merge pred and prey matrices qpress.WSS28 <- qpress.pred + qpress.prey #add/change row and column names FG<-as.vector(new.names) qpress.WSS28<-cbind(FG, qpress.WSS28) save(qpress.WSS28, file="qpress.WSS28.RData") #remove values below 10%, 20%, 30%, 40%, 50% WSS28qnm10<-qpress.WSS28 WSS28qnm10[WSS28qnm10 >=-0.10 & WSS28qnm10<= 0.10]<-0 WSS28qnm20<-qpress.WSS28 WSS28qnm20[WSS28qnm20 >=-0.20 & WSS28qnm20<= 0.20]<-0 WSS28qnm30<-qpress.WSS28 WSS28qnm30[WSS28qnm20 >=-0.30 & WSS28qnm20<= 0.30]<-0 WSS28qnm40<-qpress.WSS28 WSS28qnm40[WSS28qnm40 >=-0.40 & WSS28qnm40<= 0.40]<-0 WSS28qnm50<-qpress.WSS28 WSS28qnm50[WSS28qnm50 >=-0.50 & WSS28qnm20<= 0.50]<-0 WSS28qnm100<-qpress.WSS28 save(WSS28qnm10, file='WSS28qnm10.RData') save(WSS28qnm20, file='WSS28qnm20.RData') save(WSS28qnm30, file='WSS28qnm30.RData') save(WSS28qnm40, file='WSS28qnm40.RData') save(WSS28qnm50, file='WSS28qnm50.RData') save(WSS28qnm100, file='WSS28qnm100.RData')
Adjacency to digraph functions
MM2Qpress <- function(data){ library(data.table) mental.sheet <- as.data.table(data) names(mental.sheet)[1] <- 'Box' n <- nrow(mental.sheet) box.names <- names(mental.sheet)[which(names(mental.sheet) != 'Box')] model <- c() for(i in 1:n){ pos <- which(mental.sheet[i, .SD, .SDcol = box.names] > 0) if(length(pos) > 0){ pos.interaction <- paste(mental.sheet[i, Box], '->', names(mental.sheet[i, pos + 1, with = F]), sep = '') model <- append(model, pos.interaction) } neg <- which(mental.sheet[i, .SD, .SDcol = box.names] < 0) if(length(neg) > 0){ neg.interaction <- paste(mental.sheet[i, Box], '-*', names(mental.sheet[i, neg + 1, with = F]), sep = '') model <- append(model, neg.interaction) } } return(model) } #function to create qnm models for qpress from signed digraphs with self-limitation make.qnm<-function(modmat){ q<-MM2Qpress(modmat) qnm<-dput(q) qnm.mod<-parse.digraph(qnm) qnm.model<-enforce.limitation(qnm.mod) } #function to create qnm models for qpress from signed digraphs WITHOUT self-limitation make.qnm.nolim<-function(modmat){ q<-MM2Qpress(modmat) qnm<-dput(q) qnm.mod<-parse.digraph(qnm) qnm.model<-qnm.mod }
Qpress for WSS28
#-------------------------------------------------------------------- #####PERTURBATIONS##### #to loop over multiple models #create folder with files in directory, call list of files set.seed(1984) #load RData files if required #create list mod<-list(WSS28qnm10, WSS28qnm20, WSS28qnm30, WSS28qnm40, WSS28qnm50, WSS28qnm100) #get MM matrix into qpress format, create model q.models<-vector('list', length(mod)) for(i in seq_along(mod)){ q.models[[i]]<-make.qnm(mod[[i]]) } #write dia files # for(i in seq_along(q.models)){ # #can add text below to further specify name # write.dia(q.models[[i]], paste0("WSS", i,".dia")) # } #loop through models to create results csv file and pdf barplots for(j in seq_along(q.models)){ A <- adjacency.matrix(q.models[[j]]) #A ## Function to generate the community matrix s <- community.sampler(q.models[[j]]) ## Function to check the validation condition #press <- press.validate(GB2el, # perturb=c("Forage_Fish"=-1), # monitor=c("Commercial_Pelagic_Fishery"=1)) #Press/Perturbation with decrease forage fish and no monitoring) press <- press.validate(q.models[[j]], perturb=c("Seals"=1), monitor=F) ## Function to define the perturbation scenario impact <- press.impact(q.models[[j]],perturb=c("Seals"=1)) ## Use 1000 simulations n.sims <- 1000 results <- 0 i <- 0 while(i < n.sims) { ## Randomly choose edges to retain z <- s$select(runif(1)) ## Sample community matrix W <- s$community() ## Check press condition and stability if(!(press(W) && stable.community(W))) next ## Monitor impact post press imp <- impact(W) imp[abs(imp)<1e-07]=0 #added to remove small perturbations results <- results + outer(sign(imp),-1:1,'==') i <- i+1 rownames(results) <- levels(q.models[[j]]$From) colnames(results) <- c('-','0','+') save(results, file=paste0("WSS28","results", j,"Seals_up",".RData")) write.csv(results, paste0("WSS28","results", j,"Seals_up",".csv")) # library(RColorBrewer) # pal <- brewer.pal(n=5,"Greys")[5:1] # #opar <- par(mar=c(5,10,1,1)+0.1) # prop <- results/rowSums(results) # r <- colSums(t(prop)*(-1:1)) # filename= paste("WSS28", "plots", j,"Seals_up", ".pdf") # pdf(filename) # #adjust barplot if you want it ordered by proportion # par(mar=c(5,10,1,1)+0.1) # barplot(t(prop), # horiz=T,cex.names=0.7,cex.axis=0.8,las=2,border=T, col=pal,xlab="Proportion") if(dev.cur() > 1) dev.off() # dev.off() } } #code for ordered barplot # pal <- brewer.pal(n=5,"Greys")[5:1] # opar <- par(mar=c(5,10,1,1)+0.1) # prop <- results/rowSums(results) # r <- colSums(t(prop)*(-1:1)) # filename= paste("wss", "plot", j,"phytoplankton up", ".pdf") # pdf(filename) # barplot(t(prop[order(r),]), # horiz=T,cex.names=0.6,cex.axis=0.8,las=2,border=T,col=pal,xlab="Proportion") # par(opar)
Qpress for WSS28 no self limiting
#-------------------------------------------------------------------- #####PERTURBATIONS##### #to loop over multiple models #create folder with files in directory, call list of files set.seed(1984) #load RData files if required #create list mod<-list(WSS28qnm10, WSS28qnm20, WSS28qnm30, WSS28qnm40, WSS28qnm50, WSS28qnm100) #get MM matrix into qpress format, create model q.models<-vector('list', length(mod)) for(i in seq_along(mod)){ q.models[[i]]<-make.qnm.nolim(mod[[i]]) } #write dia files # for(i in seq_along(q.models)){ # #can add text below to further specify name # write.dia(q.models[[i]], paste0("WSS", i,".dia")) # } #loop through models to create results csv file and pdf barplots for(j in seq_along(q.models)){ A <- adjacency.matrix(q.models[[j]]) #A ## Function to generate the community matrix s <- community.sampler(q.models[[j]]) ## Function to check the validation condition #press <- press.validate(GB2el, # perturb=c("Forage_Fish"=-1), # monitor=c("Commercial_Pelagic_Fishery"=1)) #Press/Perturbation with decrease forage fish and no monitoring) press <- press.validate(q.models[[j]], perturb=c("Phytoplankton"=1), monitor=F) ## Function to define the perturbation scenario impact <- press.impact(q.models[[j]],perturb=c("Phytoplankton"=1)) ## Use 1000 simulations n.sims <- 1000 results <- 0 i <- 0 while(i < n.sims) { ## Randomly choose edges to retain z <- s$select(runif(1)) ## Sample community matrix W <- s$community() ## Check press condition and stability if(!(press(W) && stable.community(W))) next ## Monitor impact post press imp <- impact(W) imp[abs(imp)<1e-07]=0 #added to remove small perturbations results <- results + outer(sign(imp),-1:1,'==') i <- i+1 rownames(results) <- levels(q.models[[j]]$From) colnames(results) <- c('-','0','+') save(results, file=paste0("WSS28nolim","results", j,"Seals_up",".RData")) write.csv(results, paste0("WSS28nolim","results", j,"Seals_up",".csv")) # library(RColorBrewer) # pal <- brewer.pal(n=5,"Greys")[5:1] # #opar <- par(mar=c(5,10,1,1)+0.1) # prop <- results/rowSums(results) # r <- colSums(t(prop)*(-1:1)) # filename= paste("WSS28", "plots", j,"Seals_up", ".pdf") # pdf(filename) # #adjust barplot if you want it ordered by proportion # par(mar=c(5,10,1,1)+0.1) # barplot(t(prop), # horiz=T,cex.names=0.7,cex.axis=0.8,las=2,border=T, col=pal,xlab="Proportion") if(dev.cur() > 1) dev.off() # dev.off() } } #code for ordered barplot # pal <- brewer.pal(n=5,"Greys")[5:1] # opar <- par(mar=c(5,10,1,1)+0.1) # prop <- results/rowSums(results) # r <- colSums(t(prop)*(-1:1)) # filename= paste("wss", "plot", j,"phytoplankton up", ".pdf") # pdf(filename) # barplot(t(prop[order(r),]), # horiz=T,cex.names=0.6,cex.axis=0.8,las=2,border=T,col=pal,xlab="Proportion") # par(opar)
visualize results
WSS28_10_phyto_up<-get(load("WSS28results1Phytoplankton_up.RData")) WSS28_10_phyto_down<-get(load("WSS28results1Phytoplankton_dtesting.RData")) WSS28_20_phyto_up<-get(load("WSS28results2Phytoplankton_up.RData")) WSS28_20_phyto_down<-get(load("WSS28results2Phytoplankton_dtesting.RData")) WSS28_30_phyto_up<-get(load("WSS28results3Phytoplankton_up.RData")) WSS28_30_phyto_down<-get(load("WSS28results3Phytoplankton_dtesting.RData")) WSS28_40_phyto_up<-get(load("WSS28results4Phytoplankton_up.RData")) WSS28_40_phyto_down<-get(load("WSS28results4Phytoplankton_dtesting.RData")) WSS28_50_phyto_up<-get(load("WSS28results5Phytoplankton_up.RData")) WSS28_50_phyto_down<-get(load("WSS28results5Phytoplankton_dtesting.RData")) WSS28_100_phyto_up<-get(load("WSS28results6Phytoplankton_up.RData")) WSS28_100_phyto_down<-get(load("WSS28results6Phytoplankton_dtesting.RData")) WSS28_10_spelagics_up<-get(load("WSS28results1Small pelagics_up.RData")) WSS28_10_spelagics_down<-get(load("WSS28results1Small pelagics_down.RData")) WSS28_20_spelagics_up<-get(load("WSS28results2Small pelagics_up.RData")) WSS28_20_spelagics_down<-get(load("WSS28results2Small pelagics_down.RData")) WSS28_30_spelagics_up<-get(load("WSS28results3Small pelagics_up.RData")) WSS28_30_spelagics_down<-get(load("WSS28results3Small pelagics_down.RData")) WSS28_40_spelagics_up<-get(load("WSS28results4Small pelagics_up.RData")) WSS28_40_spelagics_down<-get(load("WSS28results4Small pelagics_down.RData")) WSS28_50_spelagics_up<-get(load("WSS28results5Small pelagics_up.RData")) WSS28_50_spelagics_down<-get(load("WSS28results5Small pelagics_down.RData")) WSS28_100_spelagics_up<-get(load("WSS28results6Small pelagics_up.RData")) WSS28_100_spelagics_down<-get(load("WSS28results6Small pelagics_down.RData")) WSS28_10_Seals_up<-get(load("WSS28results1Seals_up.RData")) WSS28_10_Seals_down<-get(load("WSS28results1Seals_down.RData")) WSS28_20_Seals_up<-get(load("WSS28results2Seals_up.RData")) WSS28_20_Seals_down<-get(load("WSS28results2Seals_down.RData")) WSS28_30_Seals_up<-get(load("WSS28results3Seals_up.RData")) WSS28_30_Seals_down<-get(load("WSS28results3Seals_down.RData")) WSS28_40_Seals_up<-get(load("WSS28results4Seals_up.RData")) WSS28_40_Seals_down<-get(load("WSS28results4Seals_down.RData")) WSS28_50_Seals_up<-get(load("WSS28results5Seals_up.RData")) WSS28_50_Seals_down<-get(load("WSS28results5Seals_down.RData")) WSS28_100_Seals_up<-get(load("WSS28results6Seals_up.RData")) WSS28_100_Seals_down<-get(load("WSS28results6Seals_down.RData")) WSS28_10_Fishery_up<-get(load("WSS28results1Fishery_up.RData")) WSS28_10_Fishery_down<-get(load("WSS28results1Fishery_down.RData")) WSS28_20_Fishery_up<-get(load("WSS28results2Fishery_up.RData")) WSS28_20_Fishery_down<-get(load("WSS28results2Fishery_down.RData")) WSS28_30_Fishery_up<-get(load("WSS28results3Fishery_up.RData")) WSS28_30_Fishery_down<-get(load("WSS28results3Fishery_down.RData")) WSS28_40_Fishery_up<-get(load("WSS28results4Fishery_up.RData")) WSS28_40_Fishery_down<-get(load("WSS28results4Fishery_down.RData")) WSS28_50_Fishery_up<-get(load("WSS28results5Fishery_up.RData")) WSS28_50_Fishery_down<-get(load("WSS28results5Fishery_down.RData")) WSS28_100_Fishery_up<-get(load("WSS28results6Fishery_up.RData")) WSS28_100_Fishery_down<-get(load("WSS28results6Fishery_down.RData")) WSS28.qnm.list<-list() files <- list.files(full.names = T, recursive = T, pattern = "biomass_annual.csv") subfiles<-files[c(T,F)] nfiles<-length(subfiles) #create empty dataframe with rows 1:49 (change if using a different model with differnt number of FGs) biomass_mc<-data.frame()[1:49,] for(i in 1:nfiles){ biomass<-read.csv(paste(subfiles[[i]]),skip=9,header=TRUE) #pulls data from year 1, change by adjusting row number bio_y<-biomass[2,-1] rownames(bio_y)<-(paste0("X",i)) biomass_t<-t(bio_y) biomass_mc<-cbind(biomass_mc, biomass_t) }
single Qpress for WSS28
data<-WSS28qnm50 WSS28qnm50<-MM2Qpress(data) WSS28qnm50<-dput(WSS28qnm50) WSS28qnm50<-parse.digraph(WSS28qnm50) # write.dia(WSS28qnm10Model, file="WSS28test.dia", width=8, height=2, self=T) #enforce # WSS28qnm100el<-enforce.limitation(WSS28qnm100) ## Examine unweighted adjacency matrix A <- adjacency.matrix(WSS28qnm50) # A <- adjacency.matrix(WSS28qnm50el) ## Function to generate the community matrix s <- community.sampler(WSS28qnm50) # s <- community.sampler(WSS28qnm100el) ## Function to check the validation condition #press <- press.validate(GB2el, # perturb=c("Forage_Fish"=-1), # monitor=c("Commercial_Pelagic_Fishery"=1)) #Press/Perturbation with decrease forage fish and no monitoring) press <- press.validate(WSS28qnm50, perturb=c("Phytoplankton"=1), monitor=F) ## Function to define the perturbation scenario impact <- press.impact(WSS28qnm50,perturb=c("Phytoplankton"=1)) ## Use 10000 simulations n.sims <- 1000 results <- 0 i <- 0 while(i < n.sims) { ## Randomly choose edges to retain z <- s$select(runif(1)) ## Sample community matrix W <- s$community() ## Check press condition and stability if(!(press(W) && stable.community(W))) next ## Monitor impact post press imp <- impact(W) imp[abs(imp)<1e-07]=0 #added to remove small perturbations results <- results + outer(sign(imp),-1:1,'==') i <- i+1 } ## Print results rownames(results) <- levels(WSS28qnm50$From) colnames(results) <- c('-','0','+') #results ##save results WSS28qnm50_phyto_down<-results write.csv(WSS28qnm50_phyto_down, file="WSS2850nolimit_phyto_up.csv") save(WSS28qnm50_phyto_down, file="WSS2850nolimit_phyto_up.RData")
Rpath for WSS62
## Western Scotian Shelf model #------------------------------------------------------------------------------- #Required packages library(data.table); library(Rpath); library(here) #------------------------------------------------------------------------------- #User created functions #------------------------------------------------------------------------------- #Western Scotian Shelf Groups groups <- c('Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks', 'Large pelagic', 'Cod <1', 'Cod 1-3', 'Cod 4-6', 'Cod 7+', 'S Hake <25', 'S Hake 25-31', 'S Hake 31+', 'Halibut <46', 'Halibut 46-81', 'Halibut 82+', 'Pollock <49', 'Pollock 49+', 'D piscvores <40', 'D piscvores 40+', 'L benthivores <41', 'L bentivores 41+', 'Skates <49', 'Skates 49+', 'Dogfish', 'Redfish <22', 'Redfish 22+', 'A plaice <26', 'A plaice 26+', 'Flounders <30', 'Flounders 30+', 'Haddock <3', 'Haddock 3+', 'L sculpin <25', 'L sculpin 25+', 'Herring <4', 'Herring 4+', 'Other pelagic', 'Mackerel', 'Mesopelagic', 'Small-medium benthivores', 'Squids', 'Lobster', 'Large crabs', 'Small crabs', 'Shrimps', 'Scallop', 'Bivalves', 'Other molluscs', 'Other arthropoda', 'Echinoderms', 'Sessile benthic groups', 'Worms', 'Meiofauna','Gelatinous zoop', 'Macrozoop', 'Mesozoop', 'Microzoop', 'Microflora', 'Phytoplankton', 'Discards', 'Detritus', 'Fleet1') #Identify type (0 = consumer, 1 = producer, 2 = detritus, 3 = fleet) types <- c(rep(0, 59), 1, 2, 2, 3) #Identify stanza groups stgroups <- c(rep(NA, 6), rep('Cod', 4), rep('S hake', 3), rep('Halibut', 3), rep('Pollock', 2), rep('D piscvores', 2), rep('L benthivores', 2), rep('Skates', 2), NA, rep('Redfish', 2), rep('A plaice', 2), rep('Flounders', 2), rep('Haddock', 2), rep('S sculpin', 2), rep('Herring', 2), rep(NA, 26)) #Create parameter list object WSS62.params <- create.rpath.params(groups, types, stgroups) #Fill in model parameters #Biomass, Production, consumption biomass <- c(0.407075, 0.049754, 0.04383, 0.00617, 0.0260249, 0.023573, rep(NA, 3), 0.076371, rep(NA, 2), 0.0473429, rep(NA, 2), 0.017291, NA, 0.154529, NA, 0.26641, NA, 0.071238, NA, 0.073658, 2.102007, NA, 1.523083, NA, 0.101483, NA, 0.25631, NA, 0.806943, NA, 0.073649, NA, 2.25021, 0.86577, 0.535679, 0.012515, 0.037782, 0.173893, 0.301295, 0.157303, 0.965481, 1.302563, 0.646857, 63.37292, 2, 0.853, 12.63734, 26.02519, 7.345485, 4.09129, 0.5204, 41.31, 23.79555, 5.800508, 3.678861, 33.664, 0.063293, 1, NA) pb <- c(0.071, 0.18, 0.147007, 0.25, 0.18, 0.4, 1.55, 0.39, 0.82, 0.82, 1.6, 1.3, 1.3, 0.57, 0.27, 0.27, 0.47, 0.85, 0.44, 0.62, 0.62, 0.64, 0.2, 0.26, 0.1393199, 0.27, 0.23, 0.44, 0.52, 0.53, 0.59, 0.66, 0.42, 0.6, 0.5, 0.64, 0.81, 0.74, 0.583, 0.97, 1.41, 4, 0.91, 0.654, 1.31, 3, 0.342, 0.69, 0.75, 4.4, 0.33, 0.32, 1.25, 10.8, 15.51, 3.04, NA, 82.8, 104.938, 70.639, rep(NA, 3)) qb <- c(4.94, 14.5, 7.338632, 87.6, 4.78, 4.24, rep(NA, 3), 1.326, rep(NA, 2), 3.28, rep(NA, 2), 1.61, NA, 2.94, NA, 2.5, NA, 1.92, NA, 1.9, 1.739802, NA, 2.41, NA, 2.41, NA, 3.21, NA, 2.08, NA, 3.93, NA, 2.36, 3.31, 5.08, 18.9, 4.7, 11.33333, rep(NA, 12), 62.05, 19.5, 73, rep(NA, 6)) WSS62.params$model[, Biomass := biomass] WSS62.params$model[, PB := pb] WSS62.params$model[, QB := qb] #Production to Consumption WSS62.params$model[Group %in% c('Lobster', 'Large crabs', 'Small crabs', 'Shrimps', 'Scallop', 'Bivalves', 'Other molluscs', 'Other arthropoda', 'Echinoderms', 'Sessile benthic groups', 'Worms', 'Meiofauna'), ProdCons := 0.15] WSS62.params$model[Group == 'Mesozoop', ProdCons := 0.4] WSS62.params$model[Group %in% c('Microzoop', 'Microflora'), ProdCons := 0.5] #Biomass accumulation and unassimilated production # ba <- c(0.003257, 0, 0.003702, 0, -0.001952, 0, -0.000180536, -0.0111684, -0.01359565, # -0.00381855, rep(0, 6), -0.005120684, -0.00618116, -0.00398017, -0.0133205, # -0.008473957, -0.0071238, rep(0, 5), -0.002998298, -0.00710381, 0, 0, # 0.02020396, 0.05648601, -0.003061481, -0.00368245, rep(0, 7), 0.0301295, # rep(0, 3), 0.03234285, rep(0, 13), rep(NA, 3)) #Needed to turn off biomass accumulation for stanza groups as it is already baked in ba <- c(0.003257, 0, 0.003702, 0, -0.001952, rep(0, 37), 0.0301295, rep(0, 3), 0.03234285, rep(0, 13), rep(NA, 3)) WSS62.params$model[, BioAcc := ba] WSS62.params$model[, Unassim := c(rep(0.2, 55), rep(0.4, 3), 0.2, rep(0, 3), NA)] #Detrital fate WSS62.params$model[, Detritus := c(rep(1, 60), rep(0, 3))] WSS62.params$model[, Discards := c(rep(0, 62), 1)] #Landings/Discards land <- c(rep(0, 4), 0.00209077, 0.005067847, 0, 0.024145, 0.08195, 0.010009, 0.005438, 0.023001, 0.008442, 2.9E-05, 0.00025, 0.00299, 0.004407, 0.092313, 0, 0.058946, 0, 0.004406, 0, 0.000799, 0.01325133, 0, 0.05914, 0, 0.005274, 0, 0.021106, 0.000433, 0.089896, 0.00084, 0.001181, 0.230954, 0.719899, 0.019708, 0.048208, 0, 0, 0.001637, 0.203756, 0.016668, 0, 0.000162, 0.106389, 0.010353, 0.00236, 0, 0.03026, rep(0, 11), NA) disc <- c(0, 0.002995, 0, 2.4E-05, 0.000575, 0.000469, 0.000286, rep(0, 3), 9.6E-05, 0, 0, 3.2E-05, 0.000274, 0, 0.00029, 0, 0.002805, 0, 0.000385, 0, 0.003928, 0.009166, 0.009694, 0.00027, 0, 9E-06, 0, 0.000372, 0, 0.000111, 0, 0.000367, 0.000516, 0.000211, 0, 0.00056, 0.000276, 0, 0.000236, 1.9E-05, 0.001526, 0.00725, 0.000633, 8E-06, 0, 0, 4.1E-05, 0, 0.016926, 0.002943, rep(0, 10), NA) WSS62.params$model[, Fleet1 := land] WSS62.params$model[, Fleet1.disc := disc] #Diet composition #read DC from EwE6.6 DC.groups <- groups[which(!groups %in% c('Discards', 'Detritus', 'Fleet1'))] DC.data <- as.data.table(read.csv(here('data-raw', 'WSS and BOF 11-Diet composition.csv'))) setnames(DC.data, names(DC.data), c('X', 'Group', DC.groups)) DC.data[, c('X', 'Phytoplankton') := NULL] DC.data <- DC.data[1:63, ] DC.data[, Phytoplankton := as.numeric(0)] #Merge DC matrix with diet parameter object # WSS.params$diet <- cbind(WSS.params$diet[, Group], DC.data) # #Fix column names # setnames(WSS.params$diet, 1, 'Group') # setnames(WSS.params$diet, paste0('V', 1:60), DC.groups) WSS62.params$diet <- DC.data #Set multistanza parameters #Group parameters WSS62.params$stanzas$stgroups[, VBGF_Ksp := c(0.14, 0.46, 0.14, 0.16, 0.075, 0.14, 0.13, 0.06, 0.12, 0.21, 0.27, 0.12, 0.28)] WSS62.params$stanzas$stgroups[, Wmat := c(0.04, 0.22, 0.24, 0.13, 0.02, 0.15, 0.4, 0.06, 0.06, 0.17, 0.16, 0.15, 0.33)] WSS62.params$stanzas$stgroups[, BAB := c(-0.05, 0, 0, -0.04, -0.05, -0.10, 0, 0, -0.07, 0, 0.07, -0.05, 0)] #Individual stanza parameters WSS62.params$stanzas$stindiv[, First := c(0, 12, 48, 84, 0, 25, 48, 0, 36, 72, 0, 48, 0, 47, 0, 66, 0, 109, 0, 108, 0, 57, 0, 46, 0, 36, 0, 65, 0, 48)] WSS62.params$stanzas$stindiv[, Last := c(11, 47, 83, 400, 24, 47, 400, 35, 71, 400, 47, 400, 46, 400, 65, 400, 108, 400, 107, 400, 56, 400, 45, 400, 35, 400, 64, 400, 47, 400)] WSS62.params$stanzas$stindiv[, Z := c(1.55, 0.39, 0.82, 0.82, 1.6, 1.3, 1.3, 0.57, 0.27, 0.27, 0.47, 0.85, 0.44, 0.62, 0.62, 0.64, 0.2, 0.26, 0.27, 0.23, 0.44, 0.52, 0.53, 0.59, 0.66, 0.42, 0.6, 0.5, 0.64, 0.81)] WSS62.params$stanzas$stindiv[, Leading := c(F, F, F, T, F, F, T, F, F, T, rep(c(F, T), 10))] #Need to run this function to calculate all the stanza biomass/consumptions #The calculated biomasses are different than EwE because of the BA terms...will investigate further #Removed BioAcc for stanza groups and this works close enough (rounding errors) WSS62.params <- rpath.stanzas(WSS62.params) #Ecopath WSS62 <- rpath(WSS62.params, 'Western Scotian Shelf') usethis::use_data(WSS62, overwrite = TRUE)
Make CModels WSS62
#QPress------------------------------------------------------------------------- setwd("~/JCT work/WGNARS/RpathQNM/RpathQNM/WSS62model") Group <- c( 'Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks', 'Large pelagic', 'Cod less1', 'Cod 1to3', 'Cod 4to6', 'Cod 7plus', 'S Hake less25', 'S Hake 25to31', 'S Hake 31plus', 'Halibut less46', 'Halibut 46to81', 'Halibut 82plus', 'Pollock less49', 'Pollock 49plus', 'D piscvores less40', 'D piscvores 40plus', 'L benthivores less41', 'L bentivores 41plus', 'Skates less49', 'Skates 49plus', 'Dogfish', 'Redfish less22', 'Redfish 22plus', 'A plaice less26', 'A plaice 26plus', 'Flounders less30', 'Flounders 30plus', 'Haddock less3', 'Haddock 3plus', 'L sculpin less25', 'L sculpin 25plus', 'Herring less4', 'Herring 4plus', 'Other pelagic', 'Mackerel', 'Mesopelagic', 'Smallmedium benthivores', 'Squids', 'Lobster', 'Large crabs', 'Small crabs', 'Shrimps', 'Scallop', 'Bivalves', 'Other molluscs', 'Other arthropoda', 'Echinoderms', 'Sessile benthic groups', 'Worms', 'Meiofauna','Gelatinous zoop', 'Macrozoop', 'Mesozoop', 'Microzoop', 'Microflora', 'Phytoplankton', 'Discards', 'Detritus', 'Import') groupnamescol<-c('Group', 'Whales', 'Toothed cetaceans', 'Seals', 'Sea birds', 'Sharks', 'Large pelagic', 'Cod less1', 'Cod 1to3', 'Cod 4to6', 'Cod 7plus', 'S Hake less25', 'S Hake 25to31', 'S Hake 31plus', 'Halibut less46', 'Halibut 46to81', 'Halibut 82plus', 'Pollock less49', 'Pollock 49plus', 'D piscvores less40', 'D piscvores 40plus', 'L benthivores less41', 'L bentivores 41plus', 'Skates less49', 'Skates 49plus', 'Dogfish', 'Redfish less22', 'Redfish 22plus', 'A plaice less26', 'A plaice 26plus', 'Flounders less30', 'Flounders 30plus', 'Haddock less3', 'Haddock 3plus', 'L sculpin less25', 'L sculpin 25plus', 'Herring less4', 'Herring 4plus', 'Other pelagic', 'Mackerel', 'Mesopelagic', 'Smallmedium benthivores', 'Squids', 'Lobster', 'Large crabs', 'Small crabs', 'Shrimps', 'Scallop', 'Bivalves', 'Other molluscs', 'Other arthropoda', 'Echinoderms', 'Sessile benthic groups', 'Worms', 'Meiofauna','Gelatinous zoop', 'Macrozoop', 'Mesozoop', 'Microzoop', 'Microflora', 'Phytoplankton') #Prey matrix qpress.prey <- WSS62.params$diet #change names to be compatible with Qpress setnames(qpress.prey, groupnamescol) qpress.prey<-qpress.prey[,-1] qpress.prey<-cbind(qpress.prey, Group) pred.cols <- names(qpress.prey)[which(names(qpress.prey) != 'Group')] qpress.prey <- qpress.prey[Group != 'Import'] #Ensure that the matrices are square detritus.DC <- data.table(Group = qpress.prey[, Group], Detritus = 0, Discards = 0) qpress.prey <- merge(qpress.prey, detritus.DC, by = 'Group') #line columns and rows up groupnames <- as.character(qpress.prey$Group) groupnames<-c('Group', groupnames) setcolorder(qpress.prey, groupnames) #Predator M2/F matrix x <- copy(WSS62) ngroup <- x$NUM_LIVING + x$NUM_DEAD out <- data.frame(Group = x$Group[1:ngroup], type = x$type [1:ngroup], PB = x$PB [1:ngroup]) #Calculate F mortality # totcatch <- x$Catch + x$Discards totcatch <- x$Landings + x$Discards Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))]) # setnames(Fmort, 'Fleet1', 'Fishery') # Fmort <- as.data.frame(totcatch / x$Biomass[row(as.matrix(totcatch))]) # setnames(Fmort, paste('V', 1:x$NUM_GEARS, sep = ''), # paste('F.', x$Group[(ngroup +1):x$NUM_GROUPS], sep = '')) out <- cbind(out, Fmort[1:ngroup, ]) setnames(out, 'Fmort[1:ngroup, ]', 'Fishery') #Add Fishery to prey fishery <- as.data.table(out) fishery <- fishery[, list(Group, Fishery)] qpress.prey <- merge(qpress.prey, fishery, by="Group") qpress.prey<-qpress.prey[, -1]#remove Group column fishery.0 <- as.data.table(matrix(rep(0, nrow(fishery) + 1), 1, nrow(fishery) + 1)) # setnames(fishery.0, paste('V', 1:x$NUM_Groups, sep=''), groupnames) qpress.prey <- rbindlist(list(qpress.prey, fishery.0), use.names=FALSE) # #Predator M2/F matrix # x <- copy(WSS62.b) # ngroup <- x$NUM_LIVING + x$NUM_DEAD # out <- data.frame(Group = x$Group[1:ngroup], # type = x$type [1:ngroup], # PB = x$PB [1:ngroup]) #Calculate M0 M0 <- c(x$PB[1:x$NUM_LIVING] * (1 - x$EE[1:x$NUM_LIVING]), x$EE[(x$NUM_LIVING + 1):ngroup]) out <- cbind(out, M0) #Calculate M2 bio <- x$Biomass[1:x$NUM_LIVING] BQB <- bio * x$QB[1:x$NUM_LIVING] diet <- as.data.frame(x$DC) nodetrdiet <- diet[1:x$NUM_LIVING, ] detrdiet <- diet[(x$NUM_LIVING +1):ngroup, ] newcons <- nodetrdiet * BQB[col(as.matrix(nodetrdiet))] predM <- newcons / bio[row(as.matrix(newcons))] detcons <- detrdiet * BQB[col(as.matrix(detrdiet))] predM <- rbind(predM, detcons) # setnames(predM, paste('V', 1:x$NUM_LIVING, sep = ''), # paste('M2.', x$Group[1:x$NUM_LIVING], sep = '')) qpress.pred <- as.data.table(cbind(out, predM)) M2.cols <- names(qpress.pred)[which(!names(qpress.pred) %in% c('Group', 'type', 'PB', 'M0'))] qpress.pred[, mort.sum := rowSums(.SD), .SDcols = M2.cols] qpress.pred <- qpress.pred[, .SD / mort.sum, .SDcols = M2.cols, by = Group] qpress.pred[is.na(qpress.pred)] <- 0 #Code for thresholds # for(ipred in seq_along(M2.cols)){ # setnames(qpress.pred, M2.cols[ipred], 'V1') # qpress.pred[V1 < .1, V1 := 0] # setnames(qpress.pred, 'V1', M2.cols[ipred]) # } #make qpress.pred negative qpress.pred <- qpress.pred[, .SD * -1, .SDcol = M2.cols, by = Group] #Ensure that the matrices are square detritus.m2 <- data.table(qpress.pred[, 'Group'], M2.Detritus = 0, M2.Discards = 0) qpress.pred <- merge(qpress.pred, detritus.m2, by = 'Group') qpress.pred <- qpress.pred[, .SD, .SDcols = c(M2.cols, 'M2.Detritus', 'M2.Discards')] prednames <- as.character(c('Fishery', detritus.DC$Group)) colnames(qpress.pred) <- prednames new.names<-colnames(qpress.prey) setcolorder(qpress.pred, new.names) qpress.pred <- t(qpress.pred) qpress.pred <- cbind(qpress.pred, rep(0, nrow(qpress.pred))) colnames(qpress.pred) <- c(new.names) #Merge pred and prey matrices qpress.WSS62 <- qpress.pred + qpress.prey #add/change row and column names FG<-as.vector(new.names) qpress.WSS62<-cbind(FG, qpress.WSS62) save(qpress.WSS62, file="qpress.WSS62.RData") #remove values below 10%, 20%, 30%, 40%, 50% WSS62qnm10<-qpress.WSS62 WSS62qnm10[WSS62qnm10 >=-0.10 & WSS62qnm10<= 0.10]<-0 WSS62qnm20<-qpress.WSS62 WSS62qnm20[WSS62qnm20 >=-0.20 & WSS62qnm20<= 0.20]<-0 WSS62qnm30<-qpress.WSS62 WSS62qnm30[WSS62qnm20 >=-0.30 & WSS62qnm20<= 0.30]<-0 WSS62qnm40<-qpress.WSS62 WSS62qnm40[WSS62qnm40 >=-0.40 & WSS62qnm40<= 0.40]<-0 WSS62qnm50<-qpress.WSS62 WSS62qnm50[WSS62qnm50 >=-0.50 & WSS62qnm20<= 0.50]<-0 WSS62qnm100<-qpress.WSS62 save(WSS62qnm10, file='WSS62qnm10.RData') save(WSS62qnm20, file='WSS62qnm20.RData') save(WSS62qnm30, file='WSS62qnm30.RData') save(WSS62qnm40, file='WSS62qnm40.RData') save(WSS62qnm50, file='WSS62qnm50.RData') save(WSS62qnm100, file='WSS62qnm100.RData')
Qpress for WSS62
setwd("~/JCT work/WGNARS/RpathQNM/RpathQNM/WSS62model") MM2Qpress <- function(data){ library(data.table) mental.sheet <- as.data.table(data) names(mental.sheet)[1] <- 'Box' n <- nrow(mental.sheet) box.names <- names(mental.sheet)[which(names(mental.sheet) != 'Box')] model <- c() for(i in 1:n){ pos <- which(mental.sheet[i, .SD, .SDcol = box.names] > 0) if(length(pos) > 0){ pos.interaction <- paste(mental.sheet[i, Box], '->', names(mental.sheet[i, pos + 1, with = F]), sep = '') model <- append(model, pos.interaction) } neg <- which(mental.sheet[i, .SD, .SDcol = box.names] < 0) if(length(neg) > 0){ neg.interaction <- paste(mental.sheet[i, Box], '-*', names(mental.sheet[i, neg + 1, with = F]), sep = '') model <- append(model, neg.interaction) } } return(model) } #function to create qnm models for qpress from signed digraphs make.qnm<-function(modmat){ q<-MM2Qpress(modmat) qnm<-dput(q) qnm.mod<-parse.digraph(qnm) qnm.model<-enforce.limitation(qnm.mod) } #-------------------------------------------------------------------- #####PERTURBATIONS##### #to loop over multiple models #create folder with files in directory, call list of files set.seed(1984) #get(load RData files if required #create list # mod<-list(WSS62qnm10, WSS62qnm20, WSS62qnm30, WSS62qnm40, WSS62qnm50, WSS62qnm100) mod<-list(WSS62qnm50) #get MM matrix into qpress format, create model q.models<-vector('list', length(mod)) for(i in seq_along(mod)){ q.models[[i]]<-make.qnm(mod[[i]]) } # #write dia files # for(i in seq_along(q.models)){ # #can add text below to further specify name # write.dia(q.models[[i]], paste0("WSS", i,".dia")) # } #loop through models to create results csv file and pdf barplots for(j in seq_along(q.models)){ A <- adjacency.matrix(q.models[[j]]) #A ## Function to generate the community matrix s <- community.sampler(q.models[[j]]) ## Function to check the validation condition #press <- press.validate(GB2el, # perturb=c("Forage_Fish"=-1), # monitor=c("Commercial_Pelagic_Fishery"=1)) #Press/Perturbation with decrease forage fish and no monitoring) press <- press.validate(q.models[[j]], perturb=c("Phytoplankton"=1), monitor=F) ## Function to define the perturbation scenario impact <- press.impact(q.models[[j]],perturb=c("Phytoplankton"=1)) ## Use 1000 simulations n.sims <- 1000 results <- 0 i <- 0 while(i < n.sims) { ## Randomly choose edges to retain z <- s$select(runif(1)) ## Sample community matrix W <- s$community() ## Check press condition and stability if(!(press(W) && stable.community(W))) next ## Monitor impact post press imp <- impact(W) imp[abs(imp)<1e-07]=0 #added to remove small perturbations results <- results + outer(sign(imp),-1:1,'==') i <- i+1 rownames(results) <- levels(q.models[[j]]$From) colnames(results) <- c('-','0','+') # save(results, paste0("WSS62","results", j,"Phytoplankton_up",".RData")) write.csv(results, paste0("WSS62","results", j,"Phytoplankton_up",".csv")) library(RColorBrewer) pal <- brewer.pal(n=5,"Greys")[5:1] #opar <- par(mar=c(5,10,1,1)+0.1) prop <- results/rowSums(results) r <- colSums(t(prop)*(-1:1)) filename= paste("WSS62", "plots", j,"Phytoplankton_up", ".pdf") pdf(filename) #adjust barplot if you want it ordered by proportion par(mar=c(5,10,1,1)+0.1) barplot(t(prop), horiz=T,cex.names=0.7,cex.axis=0.8,las=2,border=T, col=pal,xlab="Proportion") dev.off() } } #code for ordered barplot # pal <- brewer.pal(n=5,"Greys")[5:1] # opar <- par(mar=c(5,10,1,1)+0.1) # prop <- results/rowSums(results) # r <- colSums(t(prop)*(-1:1)) # filename= paste("wss", "plot", j,"phytoplankton up", ".pdf") # pdf(filename) # barplot(t(prop[order(r),]), # horiz=T,cex.names=0.6,cex.axis=0.8,las=2,border=T,col=pal,xlab="Proportion") # par(opar)
gather results
#gather positive simulation results WSS28.results_pos<-WSS28.results %>% filter(Outcome=="Positive") %>% unite("model_ID", Group:Direction, remove=FALSE) %>% select(!2:5)%>% pivot_longer(!model_ID, names_to = "species", values_to = "pos") #gather neutral simulation results, calculate half WSS28.results_neu<-WSS28.results %>% filter(Outcome=="Neutral")%>% unite("model_ID", Group:Direction, remove=FALSE) %>% select(!2:5)%>% pivot_longer(!model_ID, names_to = "species", values_to = "neu") %>% mutate(half_neu=neu/2) #join positive and neutral simulations WSS28.results_posneu<-left_join(WSS28.results_pos, WSS28.results_neu) %>% mutate(result=pos+half_neu) %>% # mutate(result=(positive+500)/2) %>% select(model_ID, species,result)%>% pivot_wider(names_from=species, values_from=result, names_vary='slowest') #create columns that match with QNM results WSS28.cols<-WSS28.results %>% filter(Outcome=="Positive") %>% select(1:2) WSS28.results.2<-bind_cols(WSS28.cols, WSS28.results_posneu) #align results, correct for differences in spelling WSS28.results_ecos<-WSS28.results.2 %>% mutate( Direction=case_match(Direction,"Plus" ~ "plus", "Minus" ~ "minus"), Group=case_match(Group, "Seals" ~ "Seal", .default=Group) ) %>% add_column (Model="ecosense", model_links=0) %>% select(!model_ID)%>% unite("model_ID", Group:Direction, remove=FALSE) # WSS28.results_ecos<-WSS28.results %>% # mutate( # Direction=case_match(Direction,"Plus" ~ "plus", "Minus" ~ "minus"), # Group=case_match(Group, "Seals" ~ "Seal", .default=Group) # ) %>% # select(!Outside)%>% # filter(Outcome=="Positive") %>% # add_column (Model="ecosense", model_links=0) %>% # unite("model_ID", Group:Direction, remove=FALSE) #Bind Ecosense and QNM results WSS28_results_ecosQNM<-bind_rows(WSS28.results_ecos, WSS28_QNM_results_compiled_tr) saveRDS(WSS28_results_ecosQNM, "~/JCT work/WGNARS/RpathtoQNM/RpathQNM/data/WSS28_results_ecosQNM_wfishery.RDS")
Graphing ecosense and qnm (THIS IS OUTDATED)
library(reshape2) library(ggplot2) library(gridExtra) require(tidyverse) WSS28_results_ecosQNM <- readRDS("~/JCT work/WGNARS/RpathtoQNM/RpathQNM/data/WSS28_results_ecosQNM.RDS") ecoQNMresults<-WSS28_results_ecosQNM %>% unite("scenario", Group:Direction, remove=FALSE ) %>% unite("Model_name", Model:model_links, remove=FALSE) ecoQNMresults_long<-ecoQNMresults %>% pivot_longer (cols= 8:35, names_to="node" ) %>% mutate(sign=value) %>% mutate(sign=case_when(sign >= 600 ~ "positive", sign <600 & sign >=400~ "zero", sign <400 ~ "negative")) %>% # mutate(fill=value) %>% mutate(scaled_value= (abs(value-500))) # # ecoQNMresults_long$fill[ecoQNMresults_long$fill >800] <- "strong" # ecoQNMresults_long$fill[ecoQNMresults_long$fill >600 & ecoQNMresults_long$fill <=800] <- "weak" # ecoQNMresults_long$fill[ecoQNMresults_long$fill >400 & ecoQNMresults_long$fill <=600] <- "weak" # ecoQNMresults_long$fill[ecoQNMresults_long$fill >200 & ecoQNMresults_long$fill <=400] <- "weak" # ecoQNMresults_long$fill[ecoQNMresults_long$fill <=200] <- "strong" ######################## # Plot up responses as bubble plots dat_fishplus<- ecoQNMresults_long %>% filter(scenario=="Fishery_plus") p_fishplus<- ggplot(dat_fishplus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black", "white"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black", "white")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) dat_fishminus<- ecoQNMresults_long %>% filter(scenario=="Fishery_minus") p_fishminus<- ggplot(dat_fishminus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) dat_sealminus<- ecoQNMresults_long %>% filter(scenario=="Seal_minus") p_sealminus<- ggplot(dat_sealminus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) dat_sealplus<- ecoQNMresults_long %>% filter(scenario=="Seal_plus") p_sealplus<- ggplot(dat_sealplus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) dat_smallpelagicsminus<- ecoQNMresults_long %>% filter(scenario=="Small pelagics_minus") p_smallpelagicsminus<- ggplot(dat_smallpelagicsminus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) dat_smallpelagicsplus<- ecoQNMresults_long %>% filter(scenario=="Small pelagics_plus") p_smallpelagicsplus<- ggplot(dat_smallpelagicsplus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) dat_phytoplanktonminus<- ecoQNMresults_long %>% filter(scenario=="Phytoplankton_minus") p_phytoplanktonminus<- ggplot(dat_phytoplanktonminus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) dat_phytoplanktonplus<- ecoQNMresults_long %>% filter(scenario=="Phytoplankton_plus") p_phytoplanktonplus<- ggplot(dat_phytoplanktonplus)+ geom_point(aes(y=node, x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0) ) ------------------------------------------- dat_fishminus<- ecoQNMresults_long %>% filter(scenario=="Fishery_minus") p_fishminus<- ggplot(dat_fishminus)+ geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25, 24)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + labs(title="fishery -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_sealminus<- ecoQNMresults_long %>% filter(scenario=="Seal_minus") p_sealminus<- ggplot(dat_sealminus)+ geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","black"))+ scale_shape_manual(values=c(25,4)) + scale_fill_manual(values=c("paleturquoise4", "black")) + labs(title="seal -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_sealplus<- ecoQNMresults_long %>% filter(scenario=="Seal_plus") p_sealplus<- ggplot(dat_sealplus)+ geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + labs(title="seal +")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_smallpelagicsminus<- ecoQNMresults_long %>% filter(scenario=="Small pelagics_minus") p_smallpelagicsminus<- ggplot(dat_smallpelagicsminus)+ geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + labs(title="small pelagic -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_smallpelagicsplus<- ecoQNMresults_long %>% filter(scenario=="Small pelagics_plus") p_smallpelagicsplus<- ggplot(dat_smallpelagicsplus)+ geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + labs(title="small pelagic +")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_phytoplanktonminus<- ecoQNMresults_long %>% filter(scenario=="Phytoplankton_minus") p_phytoplanktonminus<- ggplot(dat_phytoplanktonminus)+ geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + labs(title="phytoplankton -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_phytoplanktonplus<- ecoQNMresults_long %>% filter(scenario=="Phytoplankton_plus") p_phytoplanktonplus<- ggplot(dat_phytoplanktonplus)+ geom_point(aes(y=fct_inorder(node), x=Model_name,cex=scaled_value, pch=sign, color=sign, fill=sign)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c("paleturquoise4","goldenrod3","black"))+ scale_shape_manual(values=c(25,24,4)) + scale_fill_manual(values=c("paleturquoise4","goldenrod3", "black")) + labs(title="phytoplankton +")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" )
Reorganizes the combined tables with Rpath and QNM results
library(data.table); library(here) #QNM results files <- as.data.table(dir(here::here('WSS28model', 'data neutrals'))) #only use .RData files <- files[V1 %like% '.RData'] all.results <- c() for(isim in 1:nrow(files)){ load(here::here('WSS28model', 'data neutrals', files[isim, ])) #Create output data.tablereg model <- gsub(".*neu(.*)\\..*", "\\1", files[isim, ]) model_ID <- paste('QNM_', substr(model, 1, 1), 0, sep = '') scenario <- substr(model, 2, 100) groups <- dimnames(results)[[1]] neg <- results[, 1] neu <- results[, 2] pos <- results[, 3] out <- data.table(Model_ID = model_ID, Scenario = scenario, Group = groups, Neg = neg, Neu = neu, Pos = pos) #Assign Negative, Neutral, Positive, or Mixed result out[Neg > 600, Symbol := 'Neg'] out[Neu > 600, Symbol := 'Neu'] out[Pos > 600, Symbol := 'Pos'] out[is.na(Symbol), Symbol := 'Mix'] #Assign Strength of dominant direction out[Symbol == 'Neg', Size := Neg + 0.5 * Neu] out[Symbol == 'Pos', Size := Pos + 0.5 * Neu] out[Symbol == 'Neu', Size := 701] out[Symbol == 'Mix', Size := 701] out[Size > 599 & Size <= 700, Strength := 'Weak'] out[Size > 700 & Size <= 900, Strength := 'Moderate'] out[Size > 900, Strength := 'Strong'] #Drop extra columns #out[, c('Neg', 'Neu', 'Pos') := NULL] #join all.results <- rbindlist(list(all.results, out)) } #Change QNM_60 to QNM_0 all.results[Model_ID == 'QNM_60', Model_ID := 'QNM_0'] #Ecosense results load(here::here('data', 'WSS28.results.rda')) groups <- colnames(WSS28.results)[which(!colnames(WSS28.results) %in% c('Group', 'Direction', 'Outcome', 'Outside'))] #Clean up names to match QNM WSS28.results[Direction == 'Plus', Direction := 'up'] WSS28.results[Direction == 'Minus', Direction := 'down'] WSS28.results[Group == 'Small pelagics', Group := 'Spelagics'] for(isp in 1:length(groups)){ sp.result <- data.table::dcast(WSS28.results, Group + Direction ~ Outcome, value.var = groups[isp]) sp.result[, Model_ID := 'Rpath'] sp.result[, Scenario := paste(Group, Direction, sep = '_')] sp.result[, Group := NULL] sp.result[, Group := groups[isp]] #Assign Negative, Neutral, Positive, or Mixed result sp.result[Negative > 600, Symbol := 'Neg'] sp.result[Neutral > 600, Symbol := 'Neu'] sp.result[Positive > 600, Symbol := 'Pos'] sp.result[is.na(Symbol), Symbol := 'Mix'] #Assign Strength of dominant direction sp.result[Symbol == 'Neg', Size := Negative + 0.5 * Neutral] sp.result[Symbol == 'Pos', Size := Positive + 0.5 * Neutral] sp.result[Symbol == 'Neu', Size := 701] sp.result[Symbol == 'Mix', Size := 701] sp.result[Size > 599 & Size <= 700, Strength := 'Weak'] sp.result[Size > 700 & Size <= 900, Strength := 'Moderate'] sp.result[Size > 900, Strength := 'Strong'] #Drop extra columns setnames(sp.result, c('Negative', 'Neutral', 'Positive'), c('Neg', 'Neu', 'Pos')) sp.result[, 'Direction' := NULL] #join all.results <- rbindlist(list(all.results, sp.result), use.names = T) } data.table::setkey(all.results, 'Model_ID', 'Scenario', 'Group') usethis::use_data(all.results, overwrite = T)
Graphing combined plots
library(reshape2) library(ggplot2) library(gridExtra) require(tidyverse) require(grid) library(here) require(stringr) load("~/JCT work/WGNARS/RpathtoQNM/RpathQNM/data/all.results.rda") group_order<-c("Detritus", "Discards", "Phytoplankton", "Microflora", "Microzoop", "Mesozoop", "Macrozoop", "Gelatinous zoop", "Meiofauna", "Worms", "Other molluscs", "Bivalves", "Shrimps", "Small crab other arthropoda", "Megabenthos", "Squids", "Other pelagic mesopelagic", "Small pelagics", "Skates", "L benthivores", "D benthivores", "D piscivores", "Large pelagics", "Sharks", "Sea birds", "Seals", "Toothed cetaceans", "Whales", "Fishery") model_order<-c("Rpath","QNM_0", "QNM_10", "QNM_20", "QNM_30", "QNM_40", "QNM_50") all.results_noneut<-as_tibble(all.results) #change the naming of Model_ID, this has been fixed in table code so eleminate if need be. all.results_noneut<-all.results_noneut%>% mutate_if(is.character, str_replace_all, pattern="QNM_1", replacement="QNM_10") %>% mutate_if(is.character, str_replace_all, pattern="QNM_2", replacement="QNM_20") %>% mutate_if(is.character, str_replace_all, pattern="QNM_3", replacement="QNM_30") %>% mutate_if(is.character, str_replace_all, pattern="QNM_4", replacement="QNM_40") %>% mutate_if(is.character, str_replace_all, pattern="QNM_5", replacement="QNM_50") %>% mutate_if(is.character, str_replace_all, pattern="QNM_6", replacement="QNM_0") # Plot up responses as bubble plots dat_fishplus<- all.results_noneut %>% filter(Scenario=="Fishery_up") dat_fishplus<-dat_fishplus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_fishplus<- ggplot(dat_fishplus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Fishery +")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_fishminus<- all.results_noneut %>% filter(Scenario=="Fishery_down") dat_fishminus<-dat_fishminus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_fishminus<- ggplot(dat_fishminus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Fishery -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_sealplus<- all.results_noneut %>% filter(Scenario=="Seals_up") dat_sealplus<-dat_sealplus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_sealplus<- ggplot(dat_sealplus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Seal +")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_sealminus<- all.results_noneut %>% filter(Scenario=="Seals_down") dat_sealminus<-dat_sealminus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_sealminus<- ggplot(dat_sealminus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Seals -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_smallpelagicsplus<- all.results_noneut %>% filter(Scenario=="Spelagics_up") dat_smallpelagicsplus<-dat_smallpelagicsplus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_smallpelagicsplus<- ggplot(dat_smallpelagicsplus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Small pelagics +")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_smallpelagicsminus<- all.results_noneut %>% filter(Scenario=="Spelagics_down") dat_smallpelagicsminus<-dat_smallpelagicsminus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_smallpelagicsminus<- ggplot(dat_smallpelagicsminus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Small pelagics -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_phytoplanktonplus<- all.results_noneut %>% filter(Scenario=="Phytoplankton_up") dat_phytoplanktonplus<-dat_phytoplanktonplus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_phytoplanktonplus<- ggplot(dat_phytoplanktonplus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Phytoplankton +")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) dat_phytoplanktonminus<- all.results_noneut %>% filter(Scenario=="Phytoplankton_down") dat_phytoplanktonminus<-dat_phytoplanktonminus %>% arrange(match(Group, group_order))%>% arrange(match(Model_ID, model_order)) p_phytoplanktonminus<- ggplot(dat_phytoplanktonminus)+ geom_point(aes(y=fct_inorder(Group), x=fct_inorder(Model_ID),cex=Size, pch=Symbol, color=Symbol, fill=Symbol)) + # pch=sign, color=sign, fill=fill)) + scale_color_manual( values=c( "black", "paleturquoise4","black", "goldenrod3"))+ scale_shape_manual(values=c(4, 25, 20,24)) + scale_fill_manual(values=c( "black", "paleturquoise4","black", "goldenrod3")) + labs(title="Phytoplankton -")+ theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), axis.text.x=element_text (angle=-90, hjust=0), legend.position="none" ) grid.newpage() jpeg(filename="RpathQNM_combined_plots_revisedRpath_2025.jpeg", width=20, height=15, units="in", quality=100,res=600) grid.arrange(p_phytoplanktonplus, p_smallpelagicsplus, p_sealplus, p_fishplus, p_phytoplanktonminus, p_smallpelagicsminus, p_sealminus, p_fishminus, ncol = 4)
calculate outdegrees
#for perturbed nodes in all QNM models library(tidyverse); library(here) Model_ID<-c("Rpath", "QNM_0", "QNM_10", "QNM_20", "QNM_30", "QNM_40", "QNM_50") outd_qnm_10<-(WSS28qnm10[,-1]) rownames(outd_qnm_10)<-WSS28qnm10$FG outd_qnm_10[outd_qnm_10>0]<-1 outd_qnm_10[outd_qnm_10<0]<-0 outd_qnm_20<-(t(SS28qnm20)) outd_qnm_20[outd_qnm_20>0]<-1 outd_qnm_20[outd_qnm_20<0]<-0 outd_qnm_30<-WSS28qnm30 outd_qnm_30[outd_qnm_30>0]<-1 outd_qnm_30[outd_qnm_30<0]<-0 outd_qnm_40<-WSS28qnm40 outd_qnm_40[outd_qnm_40>0]<-1 outd_qnm_40[outd_qnm_40<0]<-0 outd_qnm_50<-WSS28qnm50 outd_qnm_50[outd_qnm_50>0]<-1 outd_qnm_50[outd_qnm_50<0]<-0 outd_qnm_0<-WSS28qnm100 outd_qnm_0[outd_qnm_0>0]<-1 outd_qnm_0[outd_qnm_0<0]<-0 outd_QNM_10p<-outd_qnm_10%>% filter ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum)) outd_QNM_20p<-outd_qnm_20%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum)) outd_QNM_30p<-outd_qnm_30%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum)) outd_QNM_40p<-outd_qnm_40%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum)) outd_QNM_50p<-outd_qnm_50%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum)) outd_QNM_0p<-outd_qnm_0%>% select ("Phytoplankton", "Small pelagics", "Seals", "Fishery") %>% summarize(across(Phytoplankton:Fishery, sum)) ind_rpathQNM<-rbind(outd_QNM_0p, outd_QNM_0p, outd_QNM_10p, outd_QNM_20p, outd_QNM_30p, outd_QNM_40p, outd_QNM_50p) ind_rpathQNM<-cbind(ind_rpathQNM, Model_ID) degrees<-rep(c("indegree", "outdegree", "in_out_degree"), each=7) alldeg_rpathQNM<-rbind(ind_rpathQNM, outd_rpathQNM, deg_rpathQNM) alldeg_rpathQNM<-cbind(alldeg_rpathQNM, degrees) save(alldeg_rpathQNM, file="data/alldegrees_rpathQNM.RData")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.