## RStudio Code for Equations and Figures for:
## Heimpel George E. Heimpel, Paul K. Abram, Charlotte E. Causton, 
## Sabrina L. Celis, Moshe Coll, Ian C.W. Hardy, Marc Mangel, 
## Nicholas J. Mills, Michal Segoli

## A benefit-risk analysis for biological control introductions based on the protection of native biodiversity

## Ecological Applications

## 9 January, 2024


## The BF function (benefits function)
## Equation 4 and Figure 2a

b <- seq(1, 5, 0.01)

CVben <- 0.95
CVben1 <- 0.75
CVben2 <- 0.50
CVben3 <- 0.25
CVben4 <- 0

BVF <- (b - 1)/((b -1) + b*(1 - CVben))
BVF1 <- (b - 1)/((b -1) + b*(1 - CVben1))
BVF2 <- (b - 1)/((b -1) + b*(1 - CVben2))
BVF3 <- (b - 1)/((b -1) + b*(1 - CVben3))
BVF4 <- (b - 1)/((b -1) + b*(1 - CVben4))

par(mar = c(5, 5, 3, 4)) # bottom, left, top, right
plot(b, BVF, type = "l", lwd = 3, xlim = c(1, 5), ylim = c(0,1), 
     xlab = "b, benefit to native species",
     ylab = "BF, benefits function",
     xaxp = c(1, 5, 4),
     yaxp = c(0, 1, 5),
     las = 1,  # makes the tick labels horizontal
     cex.axis = 1, cex.lab = 1)
lines(b, BVF1, lwd = 3, col = "blue")
lines(b, BVF2, lwd = 3, col = "red")
#lines(b, BVF3, lwd = 3, col = "green")
lines(b, BVF4, lwd = 3, col = "brown")

###############################################
###############################################

## The DRF function (direct risk function)
## Equation 6, Figure 2b

Pa <- 1
Ij <- seq(0, 1, 0.01)
CVnt <- 0.95
CVnt1 <- 0.75
CVnt2 <- 0.50
CVnt3 <- 0.25
CVnt4 <- 0

DRF <- (Ij/(Ij + (1 - CVnt)))
DRF1 <- (Ij/(Ij + (1 - CVnt1)))
DRF2 <- (Ij/(Ij + (1 - CVnt2)))
DRF3 <- (Ij/(Ij + (1 - CVnt3)))
DRF4 <- (Ij/(Ij + (1 - CVnt4)))

par(mar = c(5, 5, 2, 4))
plot(Ij, DRF, type = "l", lwd = 3, xlim = c(0,1), ylim = c(0,1),
     xlab = "Ij, Impact on non-target species",
     ylab = "DRF, Direct Risks Function",
     xaxp = c(0, 1, 5),
     yaxp = c(0, 1, 5),
     las = 1,  # makes the tick labels horizontal
     cex.axis = 1, cex.lab = 1)
lines(Ij, DRF1, lwd = 3, col = "green")
lines(Ij, DRF2, lwd = 3, col = "blue")
#lines(Ij, DRF3, lwd = 3, col = "brown")
lines(Ij, DRF4, lwd = 3, col = "red")

###############################################
###############################################

## The IRF function (indirect risks function)
## Equation 7, Fig. 2c

Patar <- seq(0, 1, 0.1)

sigma <- 0.5

Itar <- 1
Itar1 <- 0.67
Itar2 <- 0.33
Itar3 <- 0.25
Itar4 <- 0

IRF <- Patar*exp(-(Itar + 1 - sigma))
IRF1 <- Patar*exp(-(Itar1 + 1 - sigma))
IRF2 <- Patar*exp(-(Itar2 + 1 - sigma))
IRF3 <- Patar*exp(-(Itar3 + 1 - sigma))
IRF4 <- Patar*exp(-(Itar4 + 1 - sigma))

par(mar = c(5, 5, 2, 4))
plot(Patar, IRF, type = 'l', lwd = 3, ylim = c(0, 0.6),
     xlab = "Patar, Probability that target species is attacked",
     ylab = "IRF, Risk of indirect effects")
lines(Patar, IRF1, type = 'l', col = "blue", lwd = 3)
lines(Patar, IRF2, type = 'l', col = "red", lwd = 3)
#lines(Patar, IRF3, type = 'l', col = "brown", lwd = 5)
lines(Patar, IRF4, type = 'l', col = "green", lwd = 3)

###############################################
###############################################

## Sensitivity Analyses shown in Fig. 3
# 7 parameters: Itar, b, CVi, Paj, Ij, CVj, IRF
# vary one at a time, keeping all else at:
# b = 3, CVi = 0.75, Paj = 0.25, Ij = 0.25, CVj = 0.75, IRF = 0.15, Itar = 0.7
# Itar = 0.7
# IRF = 0.18 , from:
IRF <- 0.5*exp(-(0.5 + 1 - 0.5)); IRF
# in which - Itar = 0.5, Patar = 0.5 and sigma = 0.5

#################
# Contour plots using ggplot #
#################

library(ggplot2)
library(reshape2)

## 1. Itar x b ##

L_Itar <- 100; L_b <- 100
Itar <- seq(from = 0, to = 1, length = L_Itar); Itar
b <- seq(from = 1, to = 5, length = L_b); b
CVi <- 0.75; Paj <- 0.25; Ij <- 0.25; CVj <- 0.75; IRF <- 0.15
BRI = array(data = 0, dim = c(L_Itar, L_b)); BRI

(for(x in 1:L_Itar)
{
  BRI[x, ] = Itar[x]*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
})
## gg plot commands below
Itar_b <- melt(BRI); Itar_b
names(Itar_b) <- c('x','y','z')
fig <- ggplot(Itar_b, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 2. Itar x Ij 

L_Itar <- 100; L_Ij <- 100
Itar <- seq(from = 0, to = 1, length = L_Itar)
Ij <- seq(from = 0, to = 1, length = L_Ij)
b <- 3; CVi <- 0.75; Paj <- 0.25; CVj <- 0.75; IRF <- 0.15
BRI = array(data = 0, dim = c(L_Itar, L_Ij))
for(x in 1:L_Itar)
{
  BRI[x, ] = Itar[x]*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
Itar_Ij <- melt(BRI)
names(Itar_Ij) <- c('x','y','z')
fig <- ggplot(Itar_Ij, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 3. Itar x CVi 

L_Itar <- 100; L_CVi <- 100
Itar <- seq(from = 0, to = 1, length = L_Itar)
CVi <- seq(from = 0, to = 1, length = L_CVi)
b <- 3; Paj <- 0.25; Ij <- 0.25; CVj <- 0.75; IRF <- 0.15
BRI = array(data = 0, dim = c(L_Itar, L_CVi))
for(x in 1:L_Itar)
{
  BRI[x, ] = Itar[x]*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
Itar_CVi <- melt(BRI); Itar_CVi
names(Itar_CVi) <- c('x','y','z')
fig <- ggplot(Itar_CVi, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 4. Itar x CVj 

L_Itar <- 100; L_CVj <- 100
Itar <- seq(from = 0, to = 1, length = L_Itar)
CVj <- seq(from = 0, to = 1, length = L_CVj)
b <- 3; CVi <- 0.75; Paj <- 0.25; Ij <- 0.25; IRF <- 0.15
BRI = array(data = 0, dim = c(L_Itar, L_CVj))
for(x in 1:L_Itar)
{
  BRI[x, ] = Itar[x]*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
Itar_CVj <- melt(BRI)
names(Itar_CVj) <- c('x','y','z')
fig <- ggplot(Itar_CVj, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 5. Itar x IRF 

L_Itar <- 100; L_IRF <- 100
Itar <- seq(from = 0, to = 1, length = L_Itar)
IRF <- seq(from = 0, to = 1, length = L_IRF)
b <- 3; CVi <- 0.75; Paj <- 0.25; Ij <- 0.25; CVj <- 0.75
BRI = array(data = 0, dim = c(L_Itar, L_IRF))
for(x in 1:L_Itar)
{
  BRI[x, ] = Itar[x]*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
Itar_IRF <- melt(BRI)
names(Itar_IRF) <- c('x','y','z')
fig <- ggplot(Itar_IRF, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 6. b x Ij

L_b <- 100; L_Ij <- 100
b <- seq(from = 1, to = 5, length = L_b)
Ij <- seq(from = 0, to = 1, length = L_Ij)
Itar <- 0.7; Paj <- 0.25; CVi <- 0.75; CVj <- 0.75; IRF <- 0.15
BRI = array(data = 0, dim = c(L_b, L_Ij))
for(x in 1:L_b)
{
  BRI[x, ] = Itar*(b[x] - 1)/((b[x] -1) + b[x]*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
b_Ij <- melt(BRI)
names(b_Ij) <- c('x','y','z')
fig <- ggplot(b_Ij, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 7. b x CVi

L_b <- 100; L_CVi <- 100
b <- seq(from = 1, to = 5, length = L_b)
CVi <- seq(from = 0, to = 1, length = L_CVi)
Itar <- 0.7; Paj <- 0.25; Ij <- 0.25; CVj <- 0.75; IRF <- 0.15
BRI = array(data = 0, dim = c(L_b, L_CVi))
for(x in 1:L_b)
{
  BRI[x, ] = Itar*(b[x] - 1)/((b[x] -1) + b[x]*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
b_CVi <- melt(BRI)
names(b_CVi) <- c('x','y','z')
fig <- ggplot(b_CVi, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 8. b x CVj

L_b <- 100; L_CVj <- 100
b <- seq(from = 1, to = 5, length = L_b)
CVj <- seq(from = 0, to = 1, length = L_CVj)
Itar <- 0.7; Paj <- 0.25; CVi <- 0.75; Ij <- 0.25; IRF <- 0.15
BRI = array(data = 0, dim = c(L_b, L_CVj))
for(x in 1:L_b)
{
  BRI[x, ] = Itar*(b[x] - 1)/((b[x] -1) + b[x]*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
b_CVj <- melt(BRI)
names(b_CVj) <- c('x','y','z')
fig <- ggplot(b_CVj, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 9. b x IRF

L_b <- 100; L_IRF <- 100
b <- seq(from = 1, to = 5, length = L_b)
IRF <- seq(from = 0, to = 1, length = L_IRF)
Itar <- 0.7; Paj <- 0.25; CVi <- 0.75; Ij <- 0.25; CVj <- 0.75
BRI = array(data = 0, dim = c(L_b, L_IRF))
for(x in 1:L_b)
{
  BRI[x, ] = Itar*(b[x] - 1)/((b[x] -1) + b[x]*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
b_IRF <- melt(BRI)
names(b_IRF) <- c('x','y','z')
fig <- ggplot(b_IRF, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 10. Ij x CVi

L_Ij <- 100; L_CVi <- 100
Ij <- seq(from = 0, to = 1, length = L_Ij)
CVi <- seq(from = 0, to = 1, length = L_CVi)
Itar <- 0.7; b <- 3; CVj <- 0.75; Paj <- 0.25; IRF <- 0.15
BRI = array(data = 0, dim = c(L_Ij, L_CVi))
for(x in 1:L_Ij)
{
  BRI[x, ] = Itar*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij[x]/(Paj*Ij[x]+(1 - CVj))) + IRF)
}
## ggplot commands below#
Ij_CVi <- melt(BRI)
names(Ij_CVi) <- c('x','y','z')
fig <- ggplot(Ij_CVi, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 11. Ij x CVj

L_Ij <- 100; L_CVj <- 100
Ij <- seq(from = 0, to = 1, length = L_Ij)
CVj <- seq(from = 0, to = 1, length = L_CVj)
Itar <- 0.7; b <- 3; CVi <- 0.75; Paj <- 0.25; IRF <- 0.15
BRI = array(data = 0, dim = c(L_Ij, L_CVj))
for(x in 1:L_Ij)
{
  BRI[x, ] = Itar*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij[x]/(Paj*Ij[x]+(1 - CVj))) + IRF)
}
## ggplot commands below#
Ij_CVj <- melt(BRI)
names(Ij_CVj) <- c('x','y','z')
fig <- ggplot(Ij_CVj, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 12. Ij x IRF

L_Ij <- 100; L_IRF <- 100
Ij <- seq(from = 0, to = 1, length = L_Ij)
IRF <- seq(from = 0, to = 1, length = L_IRF)
Itar <- 0.7; b <- 3; CVi <- 0.75; Paj <- 0.25; CVj <- 0.75
BRI = array(data = 0, dim = c(L_Ij, L_IRF))
for(x in 1:L_Ij)
{
  BRI[x, ] = Itar*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij[x]/(Paj*Ij[x]+(1 - CVj))) + IRF)
}
## ggplot commands below#
Ij_IRF <- melt(BRI)
names(Ij_IRF) <- c('x','y','z')
fig <- ggplot(Ij_IRF, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 13. CVi x CVj

L_CVi <- 100; L_CVj <- 100
CVi <- seq(from = 0, to = 1, length = L_CVi)
CVj <- seq(from = 0, to = 1, length = L_CVj)
Itar <- 0.7; b <- 3; Paj <- 0.25; Ij <- 0.75; IRF <- 0.15
BRI = array(data = 0, dim = c(L_CVi, L_CVj))
for(x in 1:L_CVi)
{
  BRI[x, ] = Itar*(b - 1)/((b -1) + b*(1 - CVi[x])) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
CVi_CVj <- melt(BRI)
names(CVi_CVj) <- c('x','y','z')
fig <- ggplot(CVi_CVj, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 14. CVi x IRF

L_CVi <- 100; L_IRF <- 100
CVi <- seq(from = 0, to = 1, length = L_CVi)
IRF <- seq(from = 0, to = 1, length = L_IRF)
Itar <- 0.7; b <- 3; Paj <- 0.25; Ij <- 0.25; CVj <- 0.75
BRI = array(data = 0, dim = c(L_CVi, L_IRF))
for(x in 1:L_CVi)
{
  BRI[x, ] = Itar*(b - 1)/((b -1) + b*(1 - CVi[x])) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj))) + IRF)
}
## ggplot commands below#
CVi_IRF <- melt(BRI)
names(CVi_IRF) <- c('x','y','z')
fig <- ggplot(CVi_IRF, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

## 15. CVj x IRF

L_CVj <- 100; L_IRF <- 100
CVj <- seq(from = 0, to = 1, length = L_CVj)
IRF <- seq(from = 0, to = 1, length = L_IRF)
Itar <- 0.7; b <- 3; CVi <- 0.75; Paj <- 0.25; Ij <- 0.25
BRI = array(data = 0, dim = c(L_CVj, L_IRF))
for(x in 1:L_CVj)
{
  BRI[x, ] = Itar*(b - 1)/((b -1) + b*(1 - CVi)) - 
    (1 + IRF)*((Paj*Ij/(Paj*Ij+(1 - CVj[x]))) + IRF)
}
## ggplot commands below#
CVj_IRF <- melt(BRI)
names(CVj_IRF) <- c('x','y','z')
fig <- ggplot(CVj_IRF, aes(x,y,z = z))
fig + geom_tile(aes(fill = z)) + scale_fill_gradient2(low = "red", high = "blue")

###############################################
###############################################

## Sensitivity Analysis shown in Fig. 4
## Parameters can be varied to obtain Figs. 4a and 4b

#Effect of b (value to beneficiary species) on BRI 
#if nt and ben species are equivalend for CV

b <- seq(1, 10, 0.01) # b = 1 means no benefit to beneficiary non-target

Itar <- 1

wv <- 0.3333 # weight for vulnerability
we <- 0.3333 # weight for ecosystem services
wc <- 0.3333 # weight for cultural significance
Cv <- 0.5 # Conservation vulnerability score
Ce <- 0.5 # Conservation ecosystem service score
Cc <- 0.5 # Conservation cultural significance score
CVi <- wv*Cv + we*Ce + wc*Cc

BF <- Itar*(b - 1)/((b -1) + b*(1 - CVi)) 

Ij <- 0.50
Psj <- 1
Pej <- 1
Pej1 <- 0.5
Paj <- Psj*Pej
Paj1 <- Psj*Pej1
wv <- 0.3333 # weight for vulnerability
we <- 0.3333 # weight for ecosystem services
wc <- 0.3333 # weight for cultural significance
Cv <- 0.5 # Conversation vulnerability score
Ce <- 0.5 # Conservation ecosystem service score
Cc <- 0.5 # Conservation cultural significance score
CVj <- wv*Cv + we*Ce + wc*Cc
r <- Paj*Ij
r1 <- Paj1*Ij
DRF <- r/(r + (1 - CVj)) 
DRF1 <- r1/(r1 + (1 - CVj))
Patar <- 1 
sigma <- 0.5
IRF <- Patar*(exp(-(Itar + 1 - sigma))); IRF
BRIright <- (1 + IRF)*(DRF + IRF)
BRIright1 <- (1 + IRF)*(DRF1 + IRF)
BRI <- BF - BRIright
BRI1 <- BF - BRIright1

plot1 <- plot(b, BRI, type = "l", 
              las = 1, lwd = 3, lty = 1, ylim = c(-1, 1), xlim = c(1, 5),
              xlab = "",
              ylab = "BRI, Benefit Risk Index")
lines(b, BRI1, lwd = 3, lty = 2)
abline(h = 0, lty = 3)

###############################################
###############################################

# Fig. 5

library(Hmisc) # for minor ticks (eye roll)

BF <- 5
DRF <- seq(0, 10, 0.1)
IRF <- 0; IRF1 <- 0.5; IRF2 <- 1

BRI <- BF - (DRF + IRF)*(1 + IRF)
BRI1 <- BF - (DRF + IRF1)*(1 + IRF1)
BRI2 <- BF - (DRF + IRF2)*(1 + IRF2)

plot(DRF, BRI, type = 'l', 
     las = 1, lwd = 3, ylim = c(-10, 5),
     xlab = "",
     ylab = "")
minor.tick(nx = 2, ny = 5, tick.ratio = 0.5)
lines(DRF, BRI1, lwd = 3, lty = 2)
lines(DRF, BRI2, lwd = 3, lty = 3)
abline(h = 0, lty = 2)