tests/testthat/testthat_depletion.R

In fishR-Core-Team/FSA: Simple Fisheries Stock Assessment Methods

## Results for validation tests below ----
# fishmethods's Darter data
data(darter,package="fishmethods")
# fishmethod Leslie
fishmethods::deplet(catch=darter$catch,effort=darter$effort,method="l")
cf1fm <- l.out$results[,1:2]
ci1fm <- l.out$results[,3:4]
# FSA Leslie
ex1 <- with(darter,depletion(darter$catch,darter$effort))
cf1 <- summary(ex1)
ci1 <- confint(ex1)

# fishmethod DeLury
fishmethods::deplet(catch=darter$catch,effort=darter$effort,method="d")
cf2fm <- d.out$results[,1:2]
ci2fm <- d.out$results[,3:4]
# FSA Leslie
ex2 <- with(darter,depletion(darter$catch,darter$effort,
                             method="DeLury",Ricker.mod=TRUE))
cf2 <- summary(ex2)
ci2 <- confint(ex2) 

# DeLury's Lobster Data
data(LobsterPEI,package="FSAdata")
df <- subset(LobsterPEI,day>16)
# fishmethod Leslie
fishmethods::deplet(catch=df$catch,effort=df$effort,method="l")
cf3fm <- l.out$results[,1:2]
ci3fm <- l.out$results[,3:4]
# FSA Leslie
ex3 <- with(df,depletion(catch,effort))
cf3 <- summary(ex3)
ci3 <- confint(ex3)

# fishmethod Leslie
fishmethods::deplet(catch=df$catch,effort=df$effort,method="d")
cf4fm <- d.out$results[,1:2]
ci4fm <- d.out$results[,3:4]
# FSA DeLury
ex4 <- with(df,depletion(catch,effort,method="DeLury",Ricker.mod=TRUE))
cf4 <- summary(ex4)
ci4 <- confint(ex4)

# Fischler's Blue Crab data
data(BlueCrab,package="FSAdata")
# fishmethod Leslie
fishmethods::deplet(catch=BlueCrab$catch,effort=BlueCrab$effort,method="l")
cf5fm <- l.out$results[,1:2]
ci5fm <- l.out$results[,3:4]
# FSA Leslie
ex5 <- with(BlueCrab,depletion(catch,effort))
cf5 <- summary(ex5)
ci5 <- confint(ex5)

# fishmethod DeLury
fishmethods::deplet(catch=BlueCrab$catch,effort=BlueCrab$effort,method="d")
cf6fm <- d.out$results[,1:2]
ci6fm <- d.out$results[,3:4]
# FSA DeLury
ex6 <- with(BlueCrab,depletion(catch,effort,method="DeLury",Ricker.mod=TRUE))
cf6 <- summary(ex6)
ci6 <- confint(ex6)


# Omand's SMB data
# fishmethod Leslie
fishmethods::deplet(catch=SMBassLS$catch,effort=SMBassLS$effort,method="l")
cf7fm <- l.out$results[,1:2]
ci7fm <- l.out$results[,3:4]
# FSA Leslie
ex7 <- with(SMBassLS,depletion(catch,effort))
cf7 <- summary(ex7)
ci7 <- confint(ex7)
# FSA Leslie with Ricker mod
ex7r <- with(SMBassLS,depletion(catch,effort,Ricker.mod=TRUE))
cf7r <- summary(ex7r)
ci7r <- confint(ex7r)

# fishmethod DeLury
fishmethods::deplet(catch=SMBassLS$catch,effort=SMBassLS$effort,method="d")
cf8fm <- d.out$results[,1:2]
ci8fm <- d.out$results[,3:4]
# FSA DeLury
ex8 <- with(SMBassLS,depletion(catch,effort,method="DeLury",Ricker.mod=TRUE))
cf8 <- summary(ex8)
ci8 <- confint(ex8)


## Test Messages ----
test_that("depletion() messages",{
  ## bad formulae
  expect_error(depletion(catch~effort+day,data=SMBassLS),
               "only one RHS variable")
  expect_error(depletion(catch+day~effort,data=SMBassLS),
               "more than one variable on the LHS")
  ## wrong type
  expect_error(depletion(c(346,184,49),rep(7,3),method="Derek"),
               "should be one of")
  ## no efforts
  expect_error(depletion(c(346,184,49)),
               '"effort" is missing')
  ## bad types of data
  expect_error(depletion(c(346,184,"Derek"),rep(7,3)),
               "must be a numeric vector")
  expect_error(depletion(c(346,184,49),rep("Derek",3)),
               "must be a numeric vector")
  ## different vector sizes
  expect_error(depletion(c(346,184,49,57),rep(7,3)),
               "must be same length")
  expect_error(depletion(c(346,184,49),rep(7,4)),
               "must be same length")
  ## too few catches
  expect_error(depletion(c(346,184),rep(7,2)),
               "Must have at least 3 values")
  ## negative catches or non-non-negative efforts
  expect_error(depletion(c(346,184,-49),rep(7,3)),
               "must be non-negative")
  expect_error(depletion(c(346,184,49),c(7,3,-1)),
               "must be positive")
  expect_error(depletion(c(346,184,49),c(7,3,0)),
               "must be positive")
  ## zero catches with DeLury method
  expect_error(depletion(c(346,184,0),rep(7,3),method="DeLury"),
               "Can't have zero")
  ## Bad regressions (i.e., have non-significant or positive slopes)
  expect_warning(depletion(c(49,184,346),rep(7,3)),
                 "Estimates are suspect")
  expect_warning(depletion(c(346,144,341),rep(7,3)),
                 "Estimates are suspect")
  expect_warning(depletion(c(49,184,346),rep(7,3),method="DeLury"),
                 "Estimates are suspect")
  expect_warning(depletion(c(346,144,341),rep(7,3),method="DeLury"),
                 "Estimates are suspect")
  
  suppressWarnings(ex1 <- depletion(c(346,184,49),rep(1,3)))
  ## wrong type in methods
  expect_error(summary(ex1,parm="Derek"),
               "should be one of")
  expect_error(coef(ex1,parm="Derek"),
               "should be one of")
  expect_error(confint(ex1,parm="Derek"),
               "should be one of")
  expect_error(confint(ex1,conf.level=0),
               "must be between 0 and 1")
  expect_error(confint(ex1,conf.level=1),
               "must be between 0 and 1")
  expect_error(confint(ex1,conf.level="R"),
               "must be numeric")
})


## Test Output Types ----
test_that("depletion() output types",{
  # formula matches vector method
  l1f <- depletion(catch~effort,data=SMBassLS,method="Leslie")
  l1 <- depletion(SMBassLS$catch,SMBassLS$effort,method="Leslie")
  expect_equal(l1f,l1)
  # Just check vector method after this
  expect_equal(class(l1),"depletion")
  # coef()s
  l1A <- coef(l1)
  expect_true(is.vector(l1A))
  expect_equal(class(l1A),"numeric")
  expect_equal(length(l1A),2)
  expect_equal(names(l1A),c("No","q"))
  l1A <- coef(l1,as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),2)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("No","q"))
  l1A <- coef(l1,parm="No")
  expect_true(is.vector(l1A))
  expect_equal(class(l1A),"numeric")
  expect_equal(length(l1A),1)
  expect_equal(names(l1A),c("No"))
  l1A <- coef(l1,parm="No",as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),1)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),"No")
  l1A <- coef(l1,parm="lm")
  expect_true(is.vector(l1A))
  expect_equal(class(l1A),"numeric")
  expect_equal(length(l1A),2)
  expect_equal(names(l1A),c("(Intercept)","K"))
  l1A <- coef(l1,parm="lm",as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),2)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("Intercept","K"))
  # confint()s
  l1A <- confint(l1)
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(nrow(l1A),2)
  expect_equal(ncol(l1A),2)
  expect_equal(rownames(l1A),c("No","q"))
  expect_equal(colnames(l1A),c("95% LCI","95% UCI"))
  l1A <- confint(l1,as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),4)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("No_LCI","No_UCI","q_LCI","q_UCI"))
  l1A <- confint(l1,parm="No")
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(ncol(l1A),2)
  expect_equal(nrow(l1A),1)
  expect_equal(rownames(l1A),c("No"))
  expect_equal(colnames(l1A),c("95% LCI","95% UCI"))
  l1A <- confint(l1,parm="No",as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),2)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("No_LCI","No_UCI"))
  l1A <- confint(l1,parm="lm")
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(nrow(l1A),2)
  expect_equal(ncol(l1A),2)
  expect_equal(rownames(l1A),c("(Intercept)","K"))
  expect_equal(colnames(l1A),c("95% LCI","95% UCI"))
  l1A <- confint(l1,parm="lm",as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),4)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("(Intercept)_LCI","(Intercept)_UCI","K_LCI","K_UCI"))
  # confint()s with incl.est=TRUE
  l1A <- confint(l1,incl.est=TRUE)
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(nrow(l1A),2)
  expect_equal(ncol(l1A),3)
  expect_equal(rownames(l1A),c("No","q"))
  expect_equal(colnames(l1A),c("Est","95% LCI","95% UCI"))
  l1A <- confint(l1,as.df=TRUE,incl.est=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),6)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("No","No_LCI","No_UCI","q","q_LCI","q_UCI"))
  l1A <- confint(l1,parm="No",incl.est=TRUE)
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(ncol(l1A),3)
  expect_equal(nrow(l1A),1)
  expect_equal(rownames(l1A),c("No"))
  expect_equal(colnames(l1A),c("Est","95% LCI","95% UCI"))
  l1A <- confint(l1,parm="No",as.df=TRUE,incl.est=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),3)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("No","No_LCI","No_UCI"))
  l1A <- confint(l1,parm="lm",incl.est=TRUE)
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(nrow(l1A),2)
  expect_equal(ncol(l1A),3)
  expect_equal(rownames(l1A),c("(Intercept)","K"))
  expect_equal(colnames(l1A),c("Est","95% LCI","95% UCI"))
  l1A <- confint(l1,parm="lm",as.df=TRUE,incl.est=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),6)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("(Intercept)","(Intercept)_LCI","(Intercept)_UCI",
                            "K","K_LCI","K_UCI"))
  # summary()s
  l1A <- summary(l1)
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(nrow(l1A),2)
  expect_equal(ncol(l1A),2)
  expect_equal(rownames(l1A),c("No","q"))
  expect_equal(colnames(l1A),c("Estimate","Std. Err."))
  l1A <- summary(l1,as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),4)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("No","No_SE","q","q_SE"))
  l1A <- summary(l1,parm="No")
  expect_true(is.matrix(l1A))
  expect_equal(mode(l1A),"numeric")
  expect_equal(nrow(l1A),1)
  expect_equal(ncol(l1A),2)
  expect_equal(rownames(l1A),c("No"))
  expect_equal(colnames(l1A),c("Estimate","Std. Err."))
  l1A <- summary(l1,parm="No",as.df=TRUE)
  expect_true(is.data.frame(l1A))
  expect_equal(ncol(l1A),2)
  expect_equal(nrow(l1A),1)
  expect_equal(names(l1A),c("No","No_SE"))
  l1A <- summary(l1,parm="lm")
  expect_equal(class(l1A),"summary.lm")
  expect_equal(mode(l1A),"list")
  # anova()s
  l1A <- anova(l1)
  expect_equal(class(l1A),c("anova","data.frame"))
  expect_equal(mode(l1A),"list")
  # r-squared
  expect_true(is.numeric(rSquared(l1)))
  
  ## did not check as.df=TRUE with method="DeLury"
  d1 <- depletion(SMBassLS$catch,SMBassLS$effort,method="DeLury")
  expect_equal(class(d1),"depletion")
  # coef()s
  d1A <- coef(d1)
  expect_true(is.vector(d1A))
  expect_equal(class(d1A),"numeric")
  expect_equal(length(d1A),2)
  expect_equal(names(d1A),c("No","q"))
  d1A <- coef(d1,parm="No")
  expect_true(is.vector(d1A))
  expect_equal(class(d1A),"numeric")
  expect_equal(length(d1A),1)
  expect_equal(names(d1A),c("No"))
  d1A <- coef(d1,parm="lm")
  expect_true(is.vector(d1A))
  expect_equal(class(d1A),"numeric")
  expect_equal(length(d1A),2)
  # confint()s
  d1A <- confint(d1)
  expect_true(is.matrix(d1A))
  expect_equal(mode(d1A),"numeric")
  expect_equal(nrow(d1A),2)
  expect_equal(ncol(d1A),2)
  expect_equal(rownames(d1A),c("No","q"))
  expect_equal(colnames(d1A),c("95% LCI","95% UCI"))
  d1A <- confint(d1,parm="No")
  expect_true(is.matrix(d1A))
  expect_equal(mode(d1A),"numeric")
  expect_equal(length(d1A),2)
  expect_equal(rownames(d1A),c("No"))
  expect_equal(colnames(d1A),c("95% LCI","95% UCI"))
  d1A <- confint(d1,parm="lm")
  expect_true(is.matrix(d1A))
  expect_equal(mode(d1A),"numeric")
  expect_equal(nrow(d1A),2)
  expect_equal(ncol(d1A),2)
  expect_equal(colnames(d1A),c("95% LCI","95% UCI"))
  # summary()s
  d1A <- summary(d1)
  expect_true(is.matrix(d1A))
  expect_equal(mode(d1A),"numeric")
  expect_equal(nrow(d1A),2)
  expect_equal(ncol(d1A),2)
  expect_equal(rownames(d1A),c("No","q"))
  expect_equal(colnames(d1A),c("Estimate","Std. Err."))
  d1A <- summary(d1,parm="No")
  expect_true(is.matrix(d1A))
  expect_equal(mode(d1A),"numeric")
  expect_equal(nrow(d1A),1)
  expect_equal(ncol(d1A),2)
  expect_equal(rownames(d1A),c("No"))
  expect_equal(colnames(d1A),c("Estimate","Std. Err."))
  d1A <- summary(d1,parm="lm")
  expect_equal(class(d1A),"summary.lm")
  expect_equal(mode(d1A),"list")
  # anova()s
  d1A <- anova(d1)
  expect_equal(class(d1A),c("anova","data.frame"))
  expect_equal(mode(d1A),"list")
  # r-squared
  expect_true(is.numeric(rSquared(d1)))
})


## Validate Results ----
test_that("depletion() with 'Leslie' matches fishmethod's 'deplet' for darter data",{
  expect_equal(round(cf1[["No","Estimate"]],0),round(cf1fm[["N","Estimate"]],0))
  expect_equal(round(cf1[["No","Std. Err."]],1),round(cf1fm[["N","SE"]],1))
  expect_equal(round(cf1["q",],7),round(cf1fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci1["No",],1),round(ci1fm["N",],1))
  expect_equal(round(ci1["q",],6),round(ci1fm["q",],6))
})

test_that("depletion() with 'Leslie' matches fishmethod's 'deplet' for lobster data",{
  expect_equal(round(cf3[["No","Estimate"]],0),round(cf3fm[["N","Estimate"]],0))
  expect_equal(round(cf3[["No","Std. Err."]],1),round(cf3fm[["N","SE"]],1))
  expect_equal(round(cf3["q",],7),round(cf3fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci3["No",],1),round(ci3fm["N",],1))
  expect_equal(round(ci3["q",],5),round(ci3fm["q",],5))
})

test_that("depletion() with 'Leslie' matches example 7.1 (p. 299) in Seber (2002) for lobster data",{
  expect_equal(round(cf3["No","Estimate"],1),120.5)
  expect_equal(round(cf3["q","Estimate"],4),0.0074)
  ## way off (52%)
  #expect_equal(round(ci3["No",],0),c(77,327))
  ## off by only 0.0001 for the LCI
  #expect_equal(round(ci3["q",],4),c(0.0058,0.0090))
})

test_that("depletion() with 'Leslie' matches fishmethod's 'deplet' for blue crab data",{
  expect_equal(round(cf5[["No","Estimate"]],0),round(cf5fm[["N","Estimate"]],0))
  expect_equal(round(cf5[["No","Std. Err."]],1),round(cf5fm[["N","SE"]],1))
  expect_equal(round(cf5["q",],7),round(cf5fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci5["No",],1),round(ci5fm["N",],1))
  expect_equal(round(ci5["q",],6),round(ci5fm["q",],6))
})

test_that("depletion with 'Leslie' matches example 7.2 (p. 300) in Seber (2002) for blue crab",{
  ## Off by <0.3%
  #expect_equal(round(cf5["No","Estimate"],-2),330300)
  expect_equal(round(cf5["q","Estimate"],5),0.00056)
  ## off by ~2%
  #expect_equal(round(ci5["No",],-2),c(299600,373600))
  ## off by 0.00001 for the LCI
  #expect_equal(round(ci5["q",],5),c(0.00045,0.00067))
})

test_that("depletion() with 'Leslie' matches fishmethod's 'deplet' for SMB data",{
  expect_equal(round(cf7[["No","Estimate"]],0),round(cf7fm[["N","Estimate"]],0))
  expect_equal(round(cf7[["No","Std. Err."]],1),round(cf7fm[["N","SE"]],1))
  expect_equal(round(cf7["q",],7),round(cf7fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci7["No",],1),round(ci7fm["N",],1))
  expect_equal(round(ci7["q",],5),round(ci7fm["q",],5))
})

test_that("depletion() with 'Leslie' with Ricker.mod matches example 6.1 (p. 151) in Ricker (1975).",{
  expect_equal(round(cf7r["No","Estimate"],0),1078)
  # multiplied by 7 because Ricker fit on Ct rather than CPEt but ft was a constant 7
  expect_equal(round(cf7r["q","Estimate"]*7,5),0.10676)
  ## Way different ... Ricker used a method by DeLury rather than Seber (2002)
  #expect_equal(round(ci7r["No",],0),c(814,2507))
})


test_that("depletion() with 'DeLury' and Ricker.mod matches fishmethod's 'deplet'",{
  expect_equal(round(cf2[["No","Estimate"]],0),round(cf2fm[["N","Estimate"]],0))
  expect_equal(round(cf2[["No","Std. Err."]],1),round(cf2fm[["N","SE"]],1))
  expect_equal(round(cf2["q",],7),round(cf2fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci2["No",],1),round(ci2fm["N",],1))
  expect_equal(round(ci2["q",],6),round(ci2fm["q",],6))
})

test_that("depletion() with 'DeLury' and Ricker.mod matches fishmethod's 'deplet' for the lobster data",{
  expect_equal(round(cf4[["No","Estimate"]],0),round(cf4fm[["N","Estimate"]],0))
  expect_equal(round(cf4[["No","Std. Err."]],1),round(cf4fm[["N","SE"]],1))
  expect_equal(round(cf4["q",],7),round(cf4fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci4["No",],1),round(ci4fm["N",],1))
  expect_equal(round(ci4["q",],7),round(ci4fm["q",],7))
})

#test_that("depletion() with 'DeLury' matches DeLury (1947) for lobster data",{
## REGRESSION RESULTS WERE DIFFERENT
## off by 4%
#expect_equal(round(cf4["No","Estimate"],1),116.33)
## off by 7.7% (but that is fourth decimal) 
# expect_equal(round(cf4[[1,"q"]],7),0.0079835)
#})

test_that("depletion() with 'DeLury' and Ricker.mod matches fishmethod's 'deplet' for the blue crab data",{
  expect_equal(round(cf6[["No","Estimate"]],0),round(cf6fm[["N","Estimate"]],0))
  expect_equal(round(cf6[["No","Std. Err."]],1),round(cf6fm[["N","SE"]],1))
  expect_equal(round(cf6["q",],7),round(cf6fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci6["No",],1),round(ci6fm["N",],1))
  expect_equal(round(ci6["q",],5),round(ci6fm["q",],5))
})

test_that("depletion() with 'DeLury' and Ricker.mod matches fishmethod's 'deplet' for the SMB data",{
  expect_equal(round(cf8[["No","Estimate"]],0),round(cf8fm[["N","Estimate"]],0))
  expect_equal(round(cf8[["No","Std. Err."]],1),round(cf8fm[["N","SE"]],1))
  expect_equal(round(cf8["q",],7),round(cf8fm["q",],7),ignore_attr=TRUE)
  expect_equal(round(ci8["No",],1),round(ci8fm["N",],1))
  expect_equal(round(ci8["q",],6),round(ci8fm["q",],6))
})


test_that("depletion() with 'DeLury' with Ricker.mod matches example 6.2 (p. 154) in Ricker (1975).",{
  expect_equal(round(cf8["No","Estimate"],0),1150)
  expect_equal(round(cf8["q","Estimate"],5),0.01319)
})