demos/lutjDemo.R
In smartell/CatchMSY: Catch-MSY method for estimating MSY-based reference points
#lutj.R
#load the library
#If you need to obtain the catch MSY package use:
# devtools::install_github("smartell/CatchMSY",build.vignettes=TRUE)
devtools::install_github("merrillrudd/CatchMSY",build.vignettes=TRUE)
library(catchMSY)
.NSAMP <- 5000
.NCORE <- 2
.THEME <- theme_bw(12)
# Create a new stockID
lutjanid <- new_sID(id="Lutjanidae",dfile="./inst/extdata/lutj.dat")
# Use Beverton Holt factors to establish M,ah,
# Growth parameters taken for Taape (lutjanus kasmira)
# Reference:(Moralis-Nin & Ralston, 1989)
lutjanid$linf <- 34.0
lutjanid$vbk <- 0.29
lutjanid$la.cv <- 0.07
lutjanid$a <- 2.447128e-05 #units are lbs
lutjanid$b <- 3.154
lutjanid$winf <- lutjanid$a * lutjanid$linf^lutjanid$b
lutjanid$m <- 1.50 * lutjanid$vbk
lutjanid$ah <- 1.65 / lutjanid$m
lutjanid$gh <- 0.10 * lutjanid$ah
# Set default MSY and FMSY values based on data
lutjanid$msy <- median(lutjanid$data$catch)
lutjanid$fmsy <- 0.6 * lutjanid$m
# Set age-at-entry to the fishery.
lutjanid$sel50 <- 3
lutjanid$sel95 <- 4
# Add option to specify selectivity in terms of length.
# Set parameter sampling frame
lutjanid$dfPriorInfo$dist[1] = "lnorm"
lutjanid$dfPriorInfo$par1[1] = log(lutjanid$m)
lutjanid$dfPriorInfo$par2[1] = 0.05 * lutjanid$m
lutjanid$dfPriorInfo$dist[2] = "unif"
lutjanid$dfPriorInfo$par1[2] = 0.20 * lutjanid$m
lutjanid$dfPriorInfo$par2[2] = 1.50 * lutjanid$m
lutjanid$dfPriorInfo$dist[3] = "unif"
lutjanid$dfPriorInfo$par1[3] = quantile(lutjanid$data$catch,0.05)
lutjanid$dfPriorInfo$par2[3] = quantile(lutjanid$data$catch,0.95)
#|-----------------------------------------------------------------|#
#| CATCH ONLY METHOD |#
#|-----------------------------------------------------------------|#
M0 <- lutjanid
# remove biomass data for the Catch only Method
M0$data <- M0$data[,1:2]
# generate samples from parameter samples
M0 <- sample.sid(M0,.NSAMP)
# run model with each sample
M0 <- sir.sid(M0,ncores=.NCORE)
# Get MSY statistics
M0$msy.stats <- summary(M0$S[M0$idx,3])
xmax <- max(M0$S[,"msy"])*2
ymax <- 1200
hist(M0$S[M0$idx,"msy"], xlim=c(0, xmax), ylim=c(0,ymax), xlab="MSY", main="Compare MSY distributions")
#|------------------------------------------------------------|#
#| CATCH WITH BIOMASS |#
#|----------------------------------------------------------- |#
M1 <- lutjanid
M1 <- sample.sid(M1,.NSAMP)
M1 <- sir.sid(M1,ncores=.NCORE)
M1$msy.stats <- summary(M1$S[M1$idx,3])
par(new=TRUE)
hist(M1$S[M1$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#AA000050", xaxt="n", yaxt="n", xlab="", ylab="", main="")
## adjust prior to have larger upper bound on possible MSY
M2 <- lutjanid
M2$dfPriorInfo$par1[3] = quantile(lutjanid$data$catch,0.05)
M2$dfPriorInfo$par2[3] = 1.5*quantile(lutjanid$data$catch,0.95)
# generate random samples from parameter ranges
M2 <- sample.sid(M2,.NSAMP)
# run model with each sample
M2 <- sir.sid(M2,ncores=.NCORE)
# Get MSY statistics
M2$msy.stats <- summary(M2$S[M2$idx,3])
par(new=TRUE)
hist(M2$S[M2$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#33000050", xaxt="n", yaxt="n", xlab="", ylab="", main="")
#|----------------------------------------------------------- |#
#| SIMULATION MODEL |#
#|------------------------------------------------------------|#
S1 <- lutjanid
R1 <- catchMSYModel(S1)
ii <- which(!is.na(S1$data$biomass))
# xx <- R1$bt[ii]*exp(rnorm(length(ii),0,S1$data$biomass.lse[ii]))
# S1$data$biomass[ii] <- xx
S1$data$biomass[ii] <- rlnorm(length(ii),log(R1$bt[ii]),S1$data$biomass.lse[ii])
print(S1$msy)
S1 <- sample.sid(S1,.NSAMP)
S1 <- sir.sid(S1,ncores=.NCORE)
S1$msy.stats <- summary(S1$S[S1$idx,3])
#|-------------------------------------------------------|#
#| CATCH WITH SIZE COMP |#
#|-------------------------------------------------------|#
## simulate length comp data
Ssim <- lutjanid
Ssim$la.cv <- 0.10 ## adjust length CV to simulate messy data
Rsim <- catchMSYModel(Ssim)
## add generated length comp data over time to observed catch data
M3 <- lutjanid
lc <- t(Rsim$LF)
M3$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], lc)
# generate samples from parameter samples
M3 <- sample.sid(M3,.NSAMP)
# run model with each sample
M3 <- sir.sid(M3,ncores=.NCORE)
# Get MSY statistics
M3$msy.stats <- summary(M3$S[M3$idx,3])
par(new=TRUE)
hist(M3$S[M3$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#00AA0050", xaxt="n", yaxt="n", xlab="", ylab="", main="")
## add generated length comp data over time to observed catch data and biomass data
M4 <- lutjanid
lc <- t(Rsim$LF)
M4$data <- cbind(lutjanid$data, lc)
# generate samples from parameter samples
M4 <- sample.sid(M4,.NSAMP)
# run model with each sample
M4 <- sir.sid(M4,ncores=.NCORE)
# Get MSY statistics
M4$msy.stats <- summary(M4$S[M4$idx,3])
par(new=TRUE)
hist(M4$S[M4$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#0000AA50", xaxt="n", yaxt="n", xlab="", ylab="", main="")
## less than all years of length comp data
lc_sub <- lc
lc_sub[1:(nrow(lc)-1),] <- NA
M3v2 <- M3
M3v2$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], lc_sub)
M3v2 <- sample.sid(M3v2,.NSAMP)
M3v2 <- sir.sid(M3v2,ncores=.NCORE)
M3v2$msy.stats <- summary(M3v2$S[M3v2$idx,3])
par(new=TRUE)
hist(M3v2$S[M3v2$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#FFFF0050", xaxt="n", yaxt="n", xlab="", ylab="")
#|-------------------------------------------------------|#
#| CATCH WITH MEAN LENGTH |#
#|-------------------------------------------------------|#
## simulate length comp data
Ssim <- lutjanid
Ssim$la.cv <- 0.10 ## adjust length CV to simulate messy data
Rsim <- catchMSYModel(Ssim)
## add generated mean length to observed catch data
M5 <- lutjanid
M5$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], "meanlength"=Rsim$ML, "meanlength.se"=rep(0.6, length(Rsim$ML)))
# generate samples from parameter samples
M5 <- sample.sid(M5,.NSAMP)
# run model with each sample
M5 <- sir.sid(M5,ncores=.NCORE)
# Get MSY statistics
M5$msy.stats <- summary(M5$S[M5$idx,3])
par(new=TRUE)
hist(M5$S[M5$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#0000AA50", xaxt="n", yaxt="n", xlab="", ylab="")
## less than all years of mean length
ML_sub <- c(rep(NA, length(Rsim$ML)-10), Rsim$ML[(length(Rsim$ML)-9):length(Rsim$ML)])
ML_sub_se <- rep(NA, length(ML_sub))
ML_sub_se[which(is.na(ML_sub)==FALSE)] <- 0.6
M5v2 <- M5
M5v2$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], "meanlength"=ML_sub, "meanlength.se"=ML_sub_se)
M5v2 <- sample.sid(M5v2,.NSAMP)
M5v2 <- sir.sid(M5v2,ncores=.NCORE)
M5v2$msy.stats <- summary(M5v2$S[M5v2$idx,3])
par(new=TRUE)
hist(M5v2$S[M5v2$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#0000FF50", xaxt="n", yaxt="n", xlab="", ylab="")
#|-------------------------------------------------|#
#| GRAPHICS AND SUMMARY STATISTICS |#
#|-------------------------------------------------|#
plot(M0,.THEME)
plot(M1,.THEME)
quartz()
M1$sel50 <- 3
M1$sel95 <- 4
Q <- catchMSYModel(M1)$Q
matplot(Q,type="l")
quartz()
M1$sel50 <- 2
M1$sel95 <- 3
Q <- catchMSYModel(M1)$Q
matplot(Q,type="l")
smartell/CatchMSY documentation built on May 30, 2019, 3:07 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.
#lutj.R
#load the library
#If you need to obtain the catch MSY package use:
# devtools::install_github("smartell/CatchMSY",build.vignettes=TRUE)
devtools::install_github("merrillrudd/CatchMSY",build.vignettes=TRUE)
library(catchMSY)
.NSAMP <- 5000
.NCORE <- 2
.THEME <- theme_bw(12)
# Create a new stockID
lutjanid <- new_sID(id="Lutjanidae",dfile="./inst/extdata/lutj.dat")
# Use Beverton Holt factors to establish M,ah,
# Growth parameters taken for Taape (lutjanus kasmira)
# Reference:(Moralis-Nin & Ralston, 1989)
lutjanid$linf <- 34.0
lutjanid$vbk <- 0.29
lutjanid$la.cv <- 0.07
lutjanid$a <- 2.447128e-05 #units are lbs
lutjanid$b <- 3.154
lutjanid$winf <- lutjanid$a * lutjanid$linf^lutjanid$b
lutjanid$m <- 1.50 * lutjanid$vbk
lutjanid$ah <- 1.65 / lutjanid$m
lutjanid$gh <- 0.10 * lutjanid$ah
# Set default MSY and FMSY values based on data
lutjanid$msy <- median(lutjanid$data$catch)
lutjanid$fmsy <- 0.6 * lutjanid$m
# Set age-at-entry to the fishery.
lutjanid$sel50 <- 3
lutjanid$sel95 <- 4
# Add option to specify selectivity in terms of length.
# Set parameter sampling frame
lutjanid$dfPriorInfo$dist[1] = "lnorm"
lutjanid$dfPriorInfo$par1[1] = log(lutjanid$m)
lutjanid$dfPriorInfo$par2[1] = 0.05 * lutjanid$m
lutjanid$dfPriorInfo$dist[2] = "unif"
lutjanid$dfPriorInfo$par1[2] = 0.20 * lutjanid$m
lutjanid$dfPriorInfo$par2[2] = 1.50 * lutjanid$m
lutjanid$dfPriorInfo$dist[3] = "unif"
lutjanid$dfPriorInfo$par1[3] = quantile(lutjanid$data$catch,0.05)
lutjanid$dfPriorInfo$par2[3] = quantile(lutjanid$data$catch,0.95)
#|-----------------------------------------------------------------|#
#| CATCH ONLY METHOD |#
#|-----------------------------------------------------------------|#
M0 <- lutjanid
# remove biomass data for the Catch only Method
M0$data <- M0$data[,1:2]
# generate samples from parameter samples
M0 <- sample.sid(M0,.NSAMP)
# run model with each sample
M0 <- sir.sid(M0,ncores=.NCORE)
# Get MSY statistics
M0$msy.stats <- summary(M0$S[M0$idx,3])
xmax <- max(M0$S[,"msy"])*2
ymax <- 1200
hist(M0$S[M0$idx,"msy"], xlim=c(0, xmax), ylim=c(0,ymax), xlab="MSY", main="Compare MSY distributions")
#|------------------------------------------------------------|#
#| CATCH WITH BIOMASS |#
#|----------------------------------------------------------- |#
M1 <- lutjanid
M1 <- sample.sid(M1,.NSAMP)
M1 <- sir.sid(M1,ncores=.NCORE)
M1$msy.stats <- summary(M1$S[M1$idx,3])
par(new=TRUE)
hist(M1$S[M1$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#AA000050", xaxt="n", yaxt="n", xlab="", ylab="", main="")
## adjust prior to have larger upper bound on possible MSY
M2 <- lutjanid
M2$dfPriorInfo$par1[3] = quantile(lutjanid$data$catch,0.05)
M2$dfPriorInfo$par2[3] = 1.5*quantile(lutjanid$data$catch,0.95)
# generate random samples from parameter ranges
M2 <- sample.sid(M2,.NSAMP)
# run model with each sample
M2 <- sir.sid(M2,ncores=.NCORE)
# Get MSY statistics
M2$msy.stats <- summary(M2$S[M2$idx,3])
par(new=TRUE)
hist(M2$S[M2$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#33000050", xaxt="n", yaxt="n", xlab="", ylab="", main="")
#|----------------------------------------------------------- |#
#| SIMULATION MODEL |#
#|------------------------------------------------------------|#
S1 <- lutjanid
R1 <- catchMSYModel(S1)
ii <- which(!is.na(S1$data$biomass))
# xx <- R1$bt[ii]*exp(rnorm(length(ii),0,S1$data$biomass.lse[ii]))
# S1$data$biomass[ii] <- xx
S1$data$biomass[ii] <- rlnorm(length(ii),log(R1$bt[ii]),S1$data$biomass.lse[ii])
print(S1$msy)
S1 <- sample.sid(S1,.NSAMP)
S1 <- sir.sid(S1,ncores=.NCORE)
S1$msy.stats <- summary(S1$S[S1$idx,3])
#|-------------------------------------------------------|#
#| CATCH WITH SIZE COMP |#
#|-------------------------------------------------------|#
## simulate length comp data
Ssim <- lutjanid
Ssim$la.cv <- 0.10 ## adjust length CV to simulate messy data
Rsim <- catchMSYModel(Ssim)
## add generated length comp data over time to observed catch data
M3 <- lutjanid
lc <- t(Rsim$LF)
M3$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], lc)
# generate samples from parameter samples
M3 <- sample.sid(M3,.NSAMP)
# run model with each sample
M3 <- sir.sid(M3,ncores=.NCORE)
# Get MSY statistics
M3$msy.stats <- summary(M3$S[M3$idx,3])
par(new=TRUE)
hist(M3$S[M3$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#00AA0050", xaxt="n", yaxt="n", xlab="", ylab="", main="")
## add generated length comp data over time to observed catch data and biomass data
M4 <- lutjanid
lc <- t(Rsim$LF)
M4$data <- cbind(lutjanid$data, lc)
# generate samples from parameter samples
M4 <- sample.sid(M4,.NSAMP)
# run model with each sample
M4 <- sir.sid(M4,ncores=.NCORE)
# Get MSY statistics
M4$msy.stats <- summary(M4$S[M4$idx,3])
par(new=TRUE)
hist(M4$S[M4$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#0000AA50", xaxt="n", yaxt="n", xlab="", ylab="", main="")
## less than all years of length comp data
lc_sub <- lc
lc_sub[1:(nrow(lc)-1),] <- NA
M3v2 <- M3
M3v2$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], lc_sub)
M3v2 <- sample.sid(M3v2,.NSAMP)
M3v2 <- sir.sid(M3v2,ncores=.NCORE)
M3v2$msy.stats <- summary(M3v2$S[M3v2$idx,3])
par(new=TRUE)
hist(M3v2$S[M3v2$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#FFFF0050", xaxt="n", yaxt="n", xlab="", ylab="")
#|-------------------------------------------------------|#
#| CATCH WITH MEAN LENGTH |#
#|-------------------------------------------------------|#
## simulate length comp data
Ssim <- lutjanid
Ssim$la.cv <- 0.10 ## adjust length CV to simulate messy data
Rsim <- catchMSYModel(Ssim)
## add generated mean length to observed catch data
M5 <- lutjanid
M5$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], "meanlength"=Rsim$ML, "meanlength.se"=rep(0.6, length(Rsim$ML)))
# generate samples from parameter samples
M5 <- sample.sid(M5,.NSAMP)
# run model with each sample
M5 <- sir.sid(M5,ncores=.NCORE)
# Get MSY statistics
M5$msy.stats <- summary(M5$S[M5$idx,3])
par(new=TRUE)
hist(M5$S[M5$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#0000AA50", xaxt="n", yaxt="n", xlab="", ylab="")
## less than all years of mean length
ML_sub <- c(rep(NA, length(Rsim$ML)-10), Rsim$ML[(length(Rsim$ML)-9):length(Rsim$ML)])
ML_sub_se <- rep(NA, length(ML_sub))
ML_sub_se[which(is.na(ML_sub)==FALSE)] <- 0.6
M5v2 <- M5
M5v2$data <- cbind(lutjanid$data[,-grep("biomass", colnames(lutjanid$data))], "meanlength"=ML_sub, "meanlength.se"=ML_sub_se)
M5v2 <- sample.sid(M5v2,.NSAMP)
M5v2 <- sir.sid(M5v2,ncores=.NCORE)
M5v2$msy.stats <- summary(M5v2$S[M5v2$idx,3])
par(new=TRUE)
hist(M5v2$S[M5v2$idx,"msy"], xlim=c(0, xmax), ylim=c(0, ymax), col="#0000FF50", xaxt="n", yaxt="n", xlab="", ylab="")
#|-------------------------------------------------|#
#| GRAPHICS AND SUMMARY STATISTICS |#
#|-------------------------------------------------|#
plot(M0,.THEME)
plot(M1,.THEME)
quartz()
M1$sel50 <- 3
M1$sel95 <- 4
Q <- catchMSYModel(M1)$Q
matplot(Q,type="l")
quartz()
M1$sel50 <- 2
M1$sel95 <- 3
Q <- catchMSYModel(M1)$Q
matplot(Q,type="l")
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.