Bio-Economic Selection Toolbox for Marine Protected Areas - BESTMPA R Package

Remi Daigle

Thank you for your interest in the R package for the bioeconomic evaluation of MPA network design! If you haven’t already done so, please refer to the manuscript describing the model results.

Usage

If you would like to use this package to your own scenario, please fork this repository and modify as needed, or install the package directly from GitHub.

install.packages("devtools")
library(devtools)
install_github("remi-daigle/BESTMPA@gh-pages")

We kindly ask that you cite both the original manuscript and this repository in any resulting publication. We recognize that there are some aspects of the model the could be closer to reality and we look forward to seeing how you would improve it. To make changes, we recommend starting with the code in this vignette (or the index.Rmd)

If you would like to replicate our results, download the repository as is and run the simulations by running the code in index.Rmd. This script will contains all the user defined parameters, and the functions are all contained within the R/ directory if you would like to inspect and or modify them individually.

Please note: This package requires the installation of R packages before use (tidyverse, data.table, ggplot2 grid, Grid2Polygons, gridExtra, igraph, maptools, Matrix, raster, readr, rgdal, rgeos, sp, spdep, tidyr) >installation of the rgdal package is notoriously difficult, please see here and here for further information should you have any problems. ### Load Packages

require(BESTMPA)
require(rgdal)
require(rgeos)
require(raster)
require(igraph)
require(ggplot2)
require(tidyverse)
require(data.table)
require(gridExtra)
require(grid)
require(spdep)
require(Matrix)
require(Grid2Polygons)
require(plotrix)

Housekeeping

Before getting started, using the housekeeping function will prepare the environment and your file structure for modelling. This function will conveniently create a results folder, clear the environment, figures, delete previous results, and increase the memory limit if desired.

results_folder <- "BESTMPA_results"
housekeeping(results_folder=results_folder,env=TRUE,fig=TRUE,delete=FALSE,memorylimit=memory.limit())

Spatial Baselayer

We need to define the spatial domain of the model. ### Exclusive Economic Zone We’re going to use the Atlantic part of the Canadian Exclusive Economic Zone from the Maritime Boundaries Geodatabase

EEZ <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="eez_iho_union_v2")
plot(EEZ, col='blue')

Let’s Remove Canadian part of the Davis Strait for EEZ since it’s beyond cod’s normal habitat.

plot(EEZ, col='blue')
EEZ <- EEZ[EEZ$marregion!="Canadian part of the Davis Strait",]
plot(EEZ, col='red',add=T)

Load cod habitat and breeding sites

The BESTMPA package is able to exclude certain habitats, or include them by default or as a priority (i.e. that habitat will be selected before proceeding to others). We’re loading cod habitat and breeding grounds generated by georeferencing figure 2 from Lough (2004).

Habitats <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="cod_habitat")
plot(Habitats,col='green')
Breeding <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="cod_breeding")[2:15,] # remove first Breeding zone, outside EEZ
plot(Breeding,col='yellow',add=T)

Establish basic grid

Define cell_size in m, minimum size of MPA’s in number of cells, default projection and which areas to include in the grid’s dataframe

# cell size in m
cell_size <- 20000

# minimum size of MPA
cells=1

# default projection
proj  <- "+proj=utm +zone=20 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
EEZ <- spTransform(EEZ,proj)

# create grid
p <- initgrid(EEZ,cell_size,proj,areas=c(Habitats=Habitats,Breeding=Breeding))

plot(p)
plot(p[p$Habitats==1,],col='green',add=T)
plot(p[p$Breeding==1,],col='yellow',add=T)

Create protection scenarios

Define basic parameters for protection scenarios:

# load existing MPAS
oldMPA <- readOGR(dsn=paste0(getwd(),"/shapefiles"),layer="MPAs_mar")
oldMPA <- spTransform(oldMPA,CRS(proj))

#remove those smaller than cell size
oldMPA <- oldMPA[gArea(oldMPA,byid = TRUE)>=cell_size^2,]

oldMPA <- apply(gCovers(oldMPA,p,byid = TRUE),1,any)|apply(gOverlaps(oldMPA,p,byid = TRUE),1,any)

# number of replicates
replicates <- 1:100

# target protection level in proportion (e.g. 0.2 is 20% protection)
MPA_coverage <- 0.10

# fixed distance for setting MPA distance in km in fixed distance scenario
fixdist <- 75000

This is how the four scenarios were created. This portion of code is a little time consuming; therefore, I write the grid to file. If using the index.Rmd and want to generate new scenarios, remove the eval=FALSE option from the chunk below

#### set status-quo scenario ####
p <- addscenario(domain=p,included=oldMPA,replicates=replicates,name="MPA_SQ",cells=cells)

#### set maximum distance ####
p <- addscenario(domain=p,included=oldMPA,MPA_coverage=MPA_coverage,replicates=replicates,name="MPA_MD",cell_size=cell_size,cells=cells)

#### set fixed distance ####
p <- addscenario(domain=p,included=oldMPA,MPA_coverage=MPA_coverage,replicates=replicates,dist=fixdist,name="MPA_FX",cell_size=cell_size,cells=cells)

#### set targeted scenario ####
p <- addscenario(domain=p,included=oldMPA,priority=p$Habitats,excluded=p$Breeding,MPA_coverage=MPA_coverage,replicates=replicates,dist=fixdist,name="MPA_TR",cell_size=cell_size,cells=cells)

writeOGR(p,dsn="shapefiles",layer="master_polygon",driver="ESRI Shapefile",overwrite_layer = TRUE)
p <- readOGR(dsn="shapefiles",layer="master_polygon")

Or you can plot all of the scenarios statically

# protect_scen_colour <- c("purple","green","red","blue")
protect_scen_colour <- c("#D01C8B","#F1B6DA","#B8E186","#4DAC26")
protect_scen_names <- c("Status Quo", "Maximum Distance", "Fixed Distance" , "Targeted")
res <- 400
qual <- 100
textsize <- 10
# list scenarios for looping
scenarios <- names(p)[grep("MPA_",names(p))] # assumes all scenarios begin with "MPA_"
scenarios <- scenarios[order(as.numeric(substr(scenarios,8,nchar(scenarios))))]
if(!file.exists(paste0(results_folder,"/figures"))) dir.create(paste0(results_folder,"/figures"))

for(r in replicates){
    if(r==1){
        jpeg(paste0(results_folder,'/figures/f1_scenarios.jpg'),height=20,width=17,units="cm",res=res,qual=qual)
    } else {
        jpeg(paste0(results_folder,'/figures/fA',r-1,'_scenarios.jpg'),height=20,width=17,units="cm",res=res,qual=qual)
    }
    
    layout(matrix(c(1,2,3,4),nrow=2))
    par(mar=c(1,1,1,1))
    nam <- names(p)[grep(paste0("MPA_.._",r,"$"),names(p))]
    for(s in nam){
        plot(EEZ)
        plot(p[p[[s]]==1,],col=protect_scen_colour[nam==s],add=T,border='transparent')
        title(protect_scen_names[which(nam==s)])
    }
    dev.off()
}

Looping for time and scenarios

Define variable to be used in the loops.

Define time

# total model run time in years (e.g. 2001:2100 would be 100 years)
time <- 2021:2071
spinup <- 20 # number of years before "time" the model starts, results from spin-up years are not saved, all scenarios start as status quo
tot_time <- (min(time)-spinup):max(time)

Define fish growth and reproduction

# Von Bertalanffy growth model parameters (Knickle and Rose 2013)
# Lt = Linf * {1-exp(-k*(t-t0))}, where Lt is length (cm) at age t (year), Linf is the asymptotic length (cm), k is the VB growth coefficient (1/year), and t0 is the x intercept (year). Linf = 112.03 (95% CI = 10.46). k = 0.13 (95% CI = 0.021). t0 = 0.18).
Linf_mean <- 112.03
Linf_SD <- 10.46/1.96
k_mean <- 0.13
k_SD <- 0.021/1.96
t0 <- 0.18

# calculate new weight - Length-weight relationship (Knickle and Rose 2013) fish$weight <- l_to_w_int * fish$length^l_to_w_power
l_to_w_int  <- 0.000011
l_to_w_power <- 2.91

# minimum age at maturity logistic equation
age_mat_steepness <- 2.5
age_mat_sigmoid <- 4


# Fecundity (size dependent). 0.5 million eggs per kg of female
fecundity <- 0.5*10^6

# Maximum age considered in matrix
maxage <- 20

# calculate which breeding area is closest to each grid cell
breeding_near <- apply(gDistance(spTransform(Breeding,proj),p,byid = T),1,which.min)

# initial abundance
initial_abun <- 250*10^6

Define sources of natural mortality

# natural mortality (Swain & Chouinard 2008)
M <- rnorm(10000,0.5938,0.0517)
M <- M[M<=1&M>=0] # eliminate any possible values of M >1 or M <0

# larval mortality (Mountain et al. 2008)
lM <- rbeta(10000,1000,1.2) #larval mortality of 99.88% (range 98.98-99.99%)
#hist(lM);mean(lM);min(lM);max(lM)

# Beverton-Holt model for carrying capacity based recruitment mortality, carrying capacity is the mean of North American carrying capacities in Table 3 of Myers et al. 2001 (mean of log CC=-1.202222222 tonnes/km^2 SD=0.9199018667)

# Habitat carrying capacity, in kg of virtual fish per cell (4779 is the number of cells in the default grid). This could be substituted with "known" habitat carrying capacity.
CCs <- (10^rnorm(10000,-1.202222222,0.9199018667))*(cell_size^2)/1000
CCs <- CCs[CCs>0][1:4779] #enforce no negative CCs
# adjust CCs to be 0 outside habitat
CCs <- CCs*p$Habitats

Define dispersal

# larval dispersal kernels are assumed to be exponential, e_fold_larvae is the e folding scale in km (the distance at which there will be fewer settlers by a factor of e). We assume that scale to be sqrt of 2cm/s*90d (avg current velocity * PLD) because we assume that the current is like a random walk
e_fold_larvae <- 2/100000*60*60*24*90*1000

# adult dispersal kernels are also assumed to be exponential, e_fold_adult (in km) was calculated from data in Lawson & Rose 2000
e_fold_adult  <- sqrt(74.139*1000)
# minimium age for adult migration (minimum size is 50 cm) Lawson & Rose 2000
min_age_migration <- 6

# generate connectivity matrix
#adult
cma <- initcm(p,e_fold_adult,cell_size)
#larvae
cml <- initcm(p,e_fold_larvae,cell_size)

Define fisherman behaviour and management

Specify spatial distribution of fish_licenses for shore distance calculation in effort calculation

# can be spatial points or spatial polygons data frame (sp package)
fish_communities <- getData('GADM', country="CAN", level=1) # provinces
fish_communities <- fish_communities[fish_communities$NAME_1=="New Brunswick"|
                                         fish_communities$NAME_1=="Newfoundland and Labrador"|
                                         fish_communities$NAME_1=="Nova Scotia"|
                                         fish_communities$NAME_1=="Prince Edward Island"|
                                         fish_communities$ID_1==11,]    # this means Quebec, the accent make bad things happen when formatting
# fish_communities <- getBigPolys(gSimplify(fish_communities,tol=0.05,topologyPreserve = T))
fish_communities <-gSimplify(fish_communities,tol=0.05,topologyPreserve = T)

plot(fish_communities)

# number of licenses per region in fish_communities
fish_licenses <- c(866, 4714, 3002, 879, 963)

# fisheries mortality at Maximum Sustainable Yield (Mountain et al. 2008)
FMSY <- 0.28

# quota set to fraction of FMSY as per precautionary principle
FMSY_buffer <- 2/3

# number of years to use in biomass estimate
biomass_est_n_years <- 5

# minimium age caught by nets (minimum size is 38 cm) Feekings et al. 2013
min_age_catch <- 4

# calculate distances from shore
distance <- gDistance(spTransform(fish_communities,proj),p,byid = T)

# list scenarios for looping
scenarios <- names(p)[grep("MPA_",names(p))] # assumes all scenarios begin with "MPA_"
scenarios <- scenarios[order(as.numeric(substr(scenarios,8,nchar(scenarios))))]

Loop over scenario and time

If using the index.Rmd and want to generate new data, remove the eval=FALSE option from the chunk below

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish <- growth(fish,y,tot_time,initial_abun,p,maxage,recruits,M)
        
        #### reproduction ####
        recruits <- reproduction(fish,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        
        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish <- dispersal(fish,cm=cma,ages=min_age_migration:maxage)
        
        #### larval dispersal ####
        recruits <- larvaldispersal(recruits,cml,lM,CCs,fish)
        
        ############################ Harvesting #################################################
        quota <- estimatequota(fish,maxage,y,tot_time,FMSY,FMSY_buffer)
        
        fish <- harvest(fish,quota,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder,y=y,time=time)
        
    }
}

close(pb)

Define cost evaluation

########################################### cost evaluation ##############################################
# normal operating cost ($) per fisherman from (Department of Fisheries and Oceans, 2007, Table A.19 Mixed Fishery Fleet)
# labour (CAD $46587) and fuel ($9008) vary with distance, the remainder does not
fish_operating_cost_ratio <- (46587+9008)/105054 # (labour+fuel)/total
# default profitability to calibrate operating cost (Department of Fisheries and Oceans, 2007, Table A.19 and 5.7 Mixed Fishery Fleet)
Status_quo_profitability <- 166184/105054 #$ catch value/$ operating expense

# landed value for cod CAD/kg (http://www.dfo-mpo.gc.ca/stats/commercial/sea-maritimes-eng.htm)
fish_landed_value <- 1.24

# Social discount rates (choose 3)
SDR <- c(0.015,0.03,0.06)

textsize <- 10
# load fish catches
filenames <- list.files(results_folder,pattern="kg_")
catch_summary <- rbindlist(lapply(filenames,processkg,results_folder,distance,fish_landed_value))
catch_summary <- calcnetvalue(catch_summary,Status_quo_profitability,fish_operating_cost_ratio)
catch_summary <- SDRvalues(catch_summary,t=time,value=catch_summary$netvalue,SDR)
catch_summary$scenario <- ordered(catch_summary$scenario,levels=c("SQ","MD","FX","TR"))


filenames <- list.files(results_folder,pattern="fish_")
fish_summary <- rbindlist(lapply(filenames,processfish,results_folder))
fish_summary$scenario <- ordered(fish_summary$scenario,levels=c("SQ","MD","FX","TR"))

levels(catch_summary$scenario) <- protect_scen_names
levels(fish_summary$scenario) <- protect_scen_names

gathered <- catch_summary %>% 
    dplyr::select(scenario,replicate,year,starts_with("SDRcumsum")) %>% 
    gather(SDR,cumsumvalue,starts_with("SDRcumsum"))
table1 <- catch_summary %>% 
    group_by(scenario) %>% 
    summarize(Catch=paste0(round(mean(kg)/1000)," (\u00b1 ",round(std.error(kg)/1000),")"),
              Distance=paste0(round(mean(distanceGV)/1000)," (\u00b1 ",round(std.error(distanceGV)/1000),")"))

table1 <- cbind(table1,fish_summary %>% 
    group_by(scenario) %>% 
    summarize(Biomass=paste0(round(mean(biomass)/1000)," (\u00b1 ",round(std.error(biomass)/1000),")")) %>% 
    select(Biomass))

table1 <- cbind(table1,catch_summary %>% 
                    group_by(scenario,replicate) %>% 
                    summarize(Mort=sum(kg<2000000)/n()*100) %>% 
                    group_by(scenario) %>% 
                    summarize(Moratorium=paste0(round(mean(Mort),2)," (\u00b1 ",round(std.error(Mort)),")")) %>% 
                    select(Moratorium))
table1[,c(1,4,2,3,5)]
##           scenario      Biomass      Catch  Distance Moratorium
## 1       Status Quo  8597 (± 16) 2938 (± 8) 327 (± 1) 5.02 (± 0)
## 2 Maximum Distance  9887 (± 12) 3074 (± 6) 319 (± 0) 1.37 (± 0)
## 3   Fixed Distance 10768 (± 12) 3161 (± 6) 314 (± 0) 1.12 (± 0)
## 4         Targeted 11282 (± 12) 3224 (± 6) 316 (± 0) 0.98 (± 0)
########## Figure 3 - plot of biomass over time #################################
p1 <- ggplot(fish_summary,aes(x=year,y=biomass/1000,colour=scenario))+
    geom_smooth(cex=2)+
    theme_classic()+
    theme(legend.position="top",
          text=element_text(size=textsize))+
    labs(x="Year",y="Total Stock Biomass (t)")+
    scale_colour_manual(values=protect_scen_colour,
                        labels=protect_scen_names,name="")
ggsave(paste0(results_folder,'/figures/f3_biomass.jpg'),p1,height=8,width=17,units="cm",dpi=res)
## `geom_smooth()` using method = 'gam'
p1
## `geom_smooth()` using method = 'gam'

########## Figure 4 - plot of catch biomass over time #################################
p1 <- ggplot(catch_summary,aes(x=year,y=kg/1000,colour=scenario))+
    geom_smooth(cex=2)+
    theme_classic()+
    theme(legend.position="top",
          text=element_text(size=textsize))+
    labs(x="Year",y="Total catch (t)")+
    scale_colour_manual(values=protect_scen_colour,
                        labels=protect_scen_names,name="")
ggsave(paste0(results_folder,'/figures/f4_catch.jpg'),p1,height=8,width=17,units="cm",dpi=res)
## `geom_smooth()` using method = 'gam'
p1
## `geom_smooth()` using method = 'gam'

########## Figure 5 - distance from shore over time #################################
p1 <- ggplot(catch_summary,aes(x=year,y=distanceGV/1000,colour=scenario))+
    geom_smooth(cex=2)+
    theme_classic()+
    theme(legend.position="top",
          text=element_text(size=textsize))+
    labs(x="Year",y="Mean distance from shore (km)")+
    scale_colour_manual(values=protect_scen_colour,
                        labels=protect_scen_names,name="")
ggsave(paste0(results_folder,'/figures/f5_distance.jpg'),p1,height=8,width=17,units="cm",dpi=res)
## `geom_smooth()` using method = 'gam'
p1
## `geom_smooth()` using method = 'gam'

########## Figure 6 - Social Discount Rate #################################
jpeg(paste0(results_folder,'/figures/f6_SDR.jpg'),height=20,width=17,units="cm",res=res,qual=qual)
p1 <- ggplot(gathered,aes(x=year,y=cumsumvalue/10^6,colour=scenario,xlab="test"))+
        geom_smooth(cex=1.5)+
        theme_classic()+
        scale_colour_manual(values=protect_scen_colour,labels = protect_scen_names,name="")+
        facet_wrap(~SDR,ncol=1,scale="free_x")+
        theme(strip.background=element_blank(),
              legend.position="none",
              strip.text=element_text(face="bold"),
              text=element_text(size=textsize))+
        labs(x="Year",y=expression(paste("Cumulative Net Present Value (10"^"6"," CAD)",sep="")))

p2 <- ggplot(gathered[gathered$year==max(gathered$year),],aes(x=scenario,y=cumsumvalue/10^6,fill=scenario))+
    # geom_jitter(size=0.5)+
    # geom_violin(alpha=0.7)+
        geom_boxplot()+
        theme_classic()+
        scale_fill_manual(values=protect_scen_colour,labels = protect_scen_names,name="")+
        facet_wrap(~SDR,ncol=1,scale="free_x")+
        theme(strip.background=element_blank(),
              legend.position="none",
              strip.text=element_text(face="bold"),
              axis.text.x=element_text(angle=20,hjust=1),
              text=element_text(size=textsize))+
        labs(x="",y=expression(paste("Net Present Value (10"^"6"," CAD)",sep="")))
g <- ggplotGrob(p1 + theme(legend.position="bottom"))$grobs
## `geom_smooth()` using method = 'gam'
legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
lheight <- sum(legend$height)

grid.arrange(legend,arrangeGrob(p1,p2,ncol=2),heights=unit.c(lheight, unit(1, "npc") - lheight))
## `geom_smooth()` using method = 'gam'
dev.off()
## png 
##   2

# Sensitivity Analysis

results_folder_sens <- "BESTMPA_sensitivity"

housekeeping(results_folder=results_folder_sens,env=FALSE,fig=TRUE,delete=FALSE,memorylimit=memory.limit())

sens_factor <- 0.1
# list scenarios for looping
scenarios <- names(p)[grep("MPA_",names(p))] # assumes all scenarios begin with "MPA_"
#select only SQ and TR scenarios
scenarios <- scenarios[c(starts_with("MPA_SQ",vars=scenarios),starts_with("MPA_TR",vars=scenarios))]
scenarios <- scenarios[order(as.numeric(substr(scenarios,8,nchar(scenarios))))] # re-order numerically
# select only first 1 scenarios of each
scenarios <- scenarios[1:50]

fecundity

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity*(1-sens_factor),age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity*(1+sens_factor),age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="fecundity_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="fecundity_plus")

        
    }
}

close(pb)

initial_abun

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun*(1-sens_factor),p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun*(1+sens_factor),p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="initial_abun_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="initial_abun_plus")

        
    }
}

close(pb)

age_mat_sigmoid

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid*(1-sens_factor),l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid*(1+sens_factor),l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="age_mat_sigmoid_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="age_mat_sigmoid_plus")

        
    }
}

close(pb)

k_mean

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean*(1-sens_factor),k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean*(1+sens_factor),k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="k_mean_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="k_mean_plus")

        
    }
}

close(pb)

min_age_migration

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun*(1-sens_factor),p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun*(1+sens_factor),p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=(min_age_migration-1):maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=(min_age_migration+1):maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_migration_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_migration_plus")

        
    }
}

close(pb)

min_age_catch

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun*(1-sens_factor),p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun*(1+sens_factor),p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=(min_age_catch-1):maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_catch_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=(min_age_catch+1):maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="min_age_catch_plus")

        
    }
}

close(pb)

e_fold_adult

# generate connectivity matrix
#adult
cma_minus <- initcm(p,e_fold_adult*(1-sens_factor),cell_size)
cma_plus <- initcm(p,e_fold_adult*(1+sens_factor),cell_size)

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma_minus,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma_plus,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_adult_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_adult_plus")

        
    }
}

close(pb)

e_fold_larvae

# generate connectivity matrix
#larvae
cml_minus <- initcm(p,e_fold_larvae*(1-sens_factor),cell_size)
cml_plus <- initcm(p,e_fold_larvae*(1+sens_factor),cell_size)

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml_minus,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml_plus,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY,FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY,FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_larvae_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="e_fold_larvae_plus")

        
    }
}

close(pb)

FMSY

pb <- progressBar("Looping through scenarios and time", "scenario",0, 1, 0)
for(s in scenarios){
    for(y in tot_time){
        gc() # keeps memory from getting clogged
        setProgressBar(pb, (which(s==scenarios)-1)/length(scenarios)+(1/length(scenarios)*((y-min(tot_time))/length(tot_time))), label=paste("Calculating scenario:",s,"year:",y))
        ############################ Growth and reproduction #################################################
        fish_minus <- growth(fish_minus,y,tot_time,initial_abun,p,maxage,recruits_minus,M)
        fish_plus <- growth(fish_plus,y,tot_time,initial_abun,p,maxage,recruits_plus,M)

        #### reproduction ####
        recruits_minus <- reproduction(fish=fish_minus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)
        recruits_plus <- reproduction(fish=fish_plus,fecundity,age_mat_steepness,age_mat_sigmoid,l_to_w_int,l_to_w_power,Linf_mean,Linf_SD,k_mean,k_SD,t0,breeding_near,p)

        ############################ Dispersal #################################################
        #### adult dispersal ####
        fish_minus <- dispersal(fish=fish_minus,cm=cma,ages=min_age_migration:maxage)
        fish_plus <- dispersal(fish=fish_plus,cm=cma,ages=min_age_migration:maxage)

        #### larval dispersal ####
        recruits_minus <- larvaldispersal(recruits_minus,cml,lM,CCs,fish_minus)
        recruits_plus <- larvaldispersal(recruits_plus,cml,lM,CCs,fish_plus)

        ############################ Harvesting #################################################
        quota_minus <- estimatequota(fish=fish_minus,maxage,y,tot_time,FMSY*(1-sens_factor),FMSY_buffer)
        quota_plus <- estimatequota(fish=fish_plus,maxage,y,tot_time,FMSY*(1+sens_factor),FMSY_buffer)

        fish_minus <- harvest(fish=fish_minus,quota_minus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="FMSY_minus")
        fish_plus <- harvest(fish=fish_plus,quota_plus,ages=min_age_catch:maxage,distance,fish_licenses,mpabefore=p$MPA_SQ_1,mpaafter=p[[s]],writefish=TRUE,writefishcatch=TRUE,writekg=TRUE,results_folder_sens,y=y,time=time, suffix="FMSY_plus")

        
    }
}

close(pb)

Evaluate Sensitivity

# load fish catches
filenames <- list.files(results_folder_sens,pattern="kg_")
catch_summary_sens <- rbindlist(lapply(filenames,processkg,results_folder_sens,distance,fish_landed_value))
catch_summary_sens <- calcnetvalue(catch_summary_sens,Status_quo_profitability,fish_operating_cost_ratio)
catch_summary_sens <- SDRvalues(catch_summary_sens,t=time,value=catch_summary_sens$netvalue,SDR)

filenames <- list.files(results_folder_sens,pattern="fish_")
fish_summary_sens <- rbindlist(lapply(filenames,processfish,results_folder_sens))

gathered <- catch_summary_sens %>% 
    dplyr::select(scenario,replicate,year,starts_with("SDRcumsum")) %>% 
    gather(SDR,cumsumvalue,starts_with("SDRcumsum"))

radialmeans <- catch_summary_sens %>%
filter(year==max(tot_time)) %>%
group_by(scenario) %>%
summarize(meanSDR=mean(`SDRcumsum = 0.015`)) %>%
mutate(variable=gsub("^[^_]+_|_[^_]+$", "", scenario),
scenario=gsub("_.*_", "_", scenario)) %>%
spread(key=scenario,value=meanSDR) %>%
ungroup() %>%
as.data.frame()


catch_summary_means <- catch_summary %>% 
    filter(year==max(tot_time)) %>% 
    group_by(scenario) %>% 
    summarize(meancumsum = mean(`SDRcumsum = 0.015`)) %>% 
    spread(scenario,meancumsum)

catch_summary_means <- catch_summary_means[rep(1,nrow(radialmeans)),]
jpeg(paste0(results_folder,'/figures/f7_sensitivity.jpg'),height=17,width=17,units="cm",res=res,qual=qual)
radial.plot(t(cbind(catch_summary_means$`Status Quo`,catch_summary_means$Targeted)),
            labels=as.character(radialmeans[,1]),
            rp.type="p",
            radial.lim=c(min(radialmeans[,2:5])*1.3,max(radialmeans[,2:5])*1.2),
            poly.col="transparent",
            line.col=protect_scen_colour[c(1,4)],
            show.grid.labels=1,
            lwd=2,
            lty=1)

radial.plot(t(radialmeans[,2:5]),
            labels=as.character(radialmeans[,1]),
            rp.type="s",
            radial.lim=c(min(radialmeans[,2:5])*1.3,max(radialmeans[,2:5])*1.2),
            point.symbols=c("-","+","-","+"),
            point.col=protect_scen_colour[c(1,1,4,4)],
            show.grid.labels=1,
            lwd=2,
            cex=2.5,
            add=T)
legend(max(radialmeans[,2:5])*0.7,
       max(radialmeans[,2:5])*2.2,
       c("Status Quo","Status Quo 90%","Status Quo 110%","Targeted","Targeted 90%","Targeted 110%"),
       col=protect_scen_colour[c(1,1,1,4,4,4)],
       lty=c(1,NA,NA,1,NA,NA),
       pch=c(NA,"-","+",NA,"-","+"),
       lwd=2,
       pt.cex=2.5,
       bty="n")
    
dev.off()
## png 
##   2