Package 'fishmethods'

Description

Calculates annual survival (S) and instantaneous total mortality rates (Z) from age frequency by using linear regression (standard and weighted), Heincke, Chapman-Robson, Poisson GLM and GLMER methods.

Usage

agesurv(type=1, age=NULL, number=NULL, full=NULL, last=NULL, estimate=c("s","z"),
method=c("lr","he","cr","crcb","ripois","wlr","pois"), sign.est=3, sign.se=3, 
 glmer.control=glmerControl(optCtrl=list(maxfun=10000),optimizer="bobyqa"))

Arguments

Details

If type = 1, the individual age data are tabulated. The age data are then subsetted based on the full and last arguments. Most calculations follow descriptions in Seber(1982), pages 414-418. If only two ages are present, a warning message is generated and the catch curve method is not calculated. Plus groups are not allowed. Any NAs represent no estimates due to some issue with model fit like convergence. If age samples were collected via a survey using gears such as seines or trawl, or were subsampled from catch, the least biased estimators are the "pois" and "crcb" methods (Nelson, 2019).

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Seber, G. A. F. 1982. The Estimation of Animal Abundance and Related Parameters, Second Edition. The Blackburn Press, Caldwell, New Jersey. 654 pages.

Maceina, M. J. and P. W. Bettoli. 1998. Variation in largemouth bass recruitment in four mainstream impoundments of the Tennessee River. N. Am. J. Fish. Manage. 18: 990-1003.

Millar, R. B. 2015. A better estimator of mortality rate from age-frequency data. Can. J. Fish. Aquat. Sci. 72: 364-375.

Nelson, G. A. 2019. Bias in common catch-curve methods applied to age frequency data from fish surveys. ICES J. Mar. Sci. doi:10.1093/icesjms/fsz085.

Quinn, T. J. and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, New York. 542 pages.

Smith, M. W. and 5 others. 2012. Recommendations for catch-curve analysis. N. Am. J. Fish. Manage. 32: 956-967.

Examples

data(rockbass)
agesurv(age=rockbass$age,full=6)

Description

Calculates the survival and mortality estimators of Jensen (1996) where net hauls are treated as samples

Usage

agesurvcl(age = NULL, group = NULL, full = NULL, last = NULL)

Arguments

Details

The individual age data are tabulated and subsetted based on full and last. The calculations follow Jensen(1996). If only two ages are present, a warning message is generated.

Value

Matrix containing estimates of annual mortality (a), annual survival (S), and instantaneous total mortality (Z) and associated standard errors.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Jensen, A. L. 1996. Ratio estimation of mortality using catch curves. Fisheries Research 27: 61-67.

See Also

agesurv

Examples

data(Jensen)
agesurvcl(age=Jensen$age,group=Jensen$group,full=0)

Description

Creates an age-length key in numbers or proportions-at-age per size.

Usage

alk(age=NULL,size=NULL,binsize=NULL,type=1)

Arguments

Details

Create age-length keys with either numbers-at-age per size class. Records with missing size values are deleted prior to calculation. Missing ages are allowed.

Value

A table of size, total numbers at size, and numbers (or proportions)-at-age per size class.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Quinn, T. J. and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, New York. 542 pages

See Also

alkD alkss alkprop

Examples

## Not run: 
  data(pinfish) 
  with(pinfish,alk(age=round(age,0),size=sl,binsize=10))
 
## End(Not run)

Description

Calculates the D statistic (sqrt of accumulated variance among ages; Lai 1987) for a range of age sample sizes using data from an age-length key. Assumes a two-stage random sampling design with proportional or fixed allocation.

Usage

alkD(x, lss = NULL, minss = NULL, maxss = NULL, sampint = NULL, 
    allocate = 1)

Arguments

Details

Following Quinn and Deriso (1999:pages 308-309), the function calculates the D statistic (sqrt of accumulated variance among ages; Lai 1987) for a range of age sample sizes defined by minss, maxss, and sampint at a given length sample size lss. The size of an age sample at a desired level of D can be obtained by the comparison. See reference to Table 8.8, p. 314 in Quinn and Deriso.

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Quinn, T. J. and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, New York. 542 pages

Lai, H.L. 1987. Optimum allocation for estimating age composition using age-length keys. U.S. Fish. Bull. 85:179-185

See Also

alkss alkprop

Examples

data(alkdata)
alkD(alkdata,lss=1000,minss=25,maxss=1000,sampint=20,allocate=1)

Description

The alkdata data frame has 39 rows and 16 columns. The age-length key for Gulf of Hauraki snapper shown in Table 8.3 of Quinn and Deriso (1999)

Usage

alkdata

Format

This data frame contains the following columns:

len

length interval

nl

number measured in length interval

A3

number of fish aged in each age class 3 within each length interval

A4

number of fish aged in each age class 4 within each length interval

A5

number of fish aged in each age class 5 within each length interval

A6

number of fish aged in each age class 6 within each length interval

A7

number of fish aged in each age class 7 within each length interval

A8

number of fish aged in each age class 8 within each length interval

A9

number of fish aged in each age class 9 within each length interval

A10

number of fish aged in each age class 10 within each length interval

A11

number of fish aged in each age class 11 within each length interval

A12

number of fish aged in each age class 12 within each length interval

A13

number of fish aged in each age class 13 within each length interval

A14

number of fish aged in each age class 14 within each length interval

A15

number of fish aged in each age class 15 within each length interval

A16

number of fish aged in each age class 16 within each length interval

Source

Quinn, T. J. and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, NY. 542 p.

Description

Calculates proportions-at-age and standard errors from an age-length key assuming a two-stage random sampling design.

Usage

alkprop(x)

Arguments

Details

If individual fish from catches are sampled randomly for lengths and then are further subsampled for age structures, Quinn and Deriso (1999: pages 304-305) showed that the proportions of fish in each age class and corresponding standard errors can be calculated assuming a two-stage random sampling design. See reference to Table 8.4, page 308 in Quinn and Deriso.

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Quinn, T. J. and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, New York. 542 pages

See Also

alkD alkss

Examples

data(alkdata) 
alkprop(alkdata)

Description

Calculates sample sizes for age subsampling assuming a two-stage random sampling design with proportional or fixed allocation.

Usage

alkss(x, lss = NULL, cv = NULL, allocate = 1)

Arguments

Details

If individual fish from catches are sampled randomly for lengths and then are further subsampled for age structures, Quinn and Deriso (1999: pages 306-309) showed that sample sizes for age structures can be determined for proportional (the number of fish aged is selected proportional to the length frequencies) and fixed (a constant number are aged per length class) allocation assuming a two-stage random sampling design. Sample sizes are determined based on the length frequency sample size, a specified coefficient of variation, and proportional or fixed allocation. The number of age classes is calculated internally. See reference to Table 8.6, p. 312 in Quinn and Deriso.

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Quinn, T. J. and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, New York. 542 pages

See Also

alkD alkprop

Examples

data(alkdata) 
alkss(alkdata,lss=1000,cv=0.25,allocate=1)

Description

This function calculates the solar zenith, azimuth and declination angles, time at sunrise, local noon and sunset, day length, and PAR (photosynthetically available radiation, 400-700 nm) under clear skies and average atmospheric conditions (marine or continental) anywhere on the surface of the earth based on date, time, and location.

Usage

astrocalc4r(day, month, year, hour, timezone, lat, lon, 
withinput = FALSE, 
seaorland = "maritime", acknowledgment = FALSE)

Arguments

Details

Astronomical definitions are based on definitions in Meeus (2009) and Seidelmann (2006). The solar zenith angle is measured between a line drawn "straight up" from the center of the earth through the observer and a line drawn from the observer to the center of the solar disk. The zenith angle reaches its lowest daily value at local noon when the sun is highest. It reaches its maximum value at night after the sun drops below the horizon. The zenith angle and all of the solar variables calculated by astrocalc4r depend on latitude, longitude, date and time of day. For example, solar zenith angles measured at the same time of day and two different locations would differ due to differences in location. Zenith angles at the same location and two different dates or times of day also differ.

Local noon is the time of day when the sun reaches its maximum elevation and minimum solar zenith angle at the observers location. This angle occurs when the leading edge of the sun first appears above, or the trailing edge disappears below the horizon (0.83o accounts for the radius of the sun when seen from the earth and for refraction by the atmosphere). Day length is the time in hours between sunrise and sunset. Solar declination and azimuth angles describe the exact position of the sun in the sky relative to an observer based on an equatorial coordinate system (Meeus 2009). Solar declination is the angular displacement of the sun above the equatorial plane. The equation of time accounts for the relative position of the observer within the time zone and is provided because it is complicated to calculate. PAR isirradiance in lux (lx, approximately W m-2) at the surface of the earth under clear skies calculated based on the solar zenith angle and assumptions about marine or terrestrial atmospheric properties. astrocalc4r calculates PAR for wavelengths between 400-700 nm. Calculations for other wavelengths can be carried out by modifying the code to use parameters from Frouin et al. (1989). Following Frouin et al. (1989), PAR is assumed to be zero at solar zenith angles >= 90o although some sunlight may be visible in the sky when the solar zenith angle is < 108o. Angles in astrocalc4r output are in degrees although radians are used internally for calculations. Time data and results are in decimal hours (e.g. 11:30 pm = 23.5 h) local time but internal calculations are in Greenwich Mean Time (GMT). The user must specify the local time zone in terms of +/- hours relative to GMT to link local time and GMT. For example, the difference between Eastern Standard Time and GMT is -5 hours. The user must ensure that any adjustments for daylight savings time are included in the timezone value. For example, timezone=-6 for Eastern daylight time.

Value

Time of solar noon, sunrise and sunset, angles of azimuth and zenith, eqtime, declination of sun, daylight length (hours) and PAR.

Author(s)

Larry Jacobson, Alan Seaver, and Jiashen Tang NOAA National Marine Fisheries Service Northeast Fisheries Science Center, 166 Water St., Woods Hole, MA 02543

[email protected]

References

Frouin, R., Lingner, D., Gautier, C., Baker, K. and Smith, R. 1989. A simple analytical formula to compute total and photosynthetically available solar irradiance at the ocean surface under clear skies. J. Geophys. Res. 94: 9731-9742.

L. D. Jacobson, L. C. Hendrickson, and J. Tang. 2015. Solar zenith angles for biological research and an expected catch model for diel vertical migration patterns that affect stock size estimates for longfin inshore squid (Doryteuthis pealeii). Canadian Journal of Fisheries and Aquatic Sciences 72: 1329-1338.

Meeus, J. 2009. Astronomical Algorithms, 2nd Edition. Willmann-Bell, Inc., Richmond, VA. Seidelmann, P.K. 2006. Explanatory Supplement to the Astronomical Almanac. University Science Books, Sausalito, CA.

Seidelmann, P.K. 2006. Explanatory Supplement to the Astronomical Almanac. University Science Books, Sausalito, CA. This function is an R implementation of:

Jacobson L, Seaver A, Tang J. 2011. AstroCalc4R: software to calculate solar zenith angle; time at sunrise, local noon and sunset; and photosynthetically available radiation based on date, time and location. US Dept Commer, Northeast Fish Sci Cent Ref Doc. 11-14; 10 p.

Examples

astrocalc4r(day=12,month=9,year=2000,hour=12,timezone=-5,lat=40.9,lon=-110)

Description

The AtkaMack data frame has 20 rows and 4 columns. Mean length-at-age data for male and female Atka Mackerel as listed in Table 3 of Kimura (1990)

Usage

AtkaMack

Format

This data frame contains the following columns:

age

fish age

len

mean length of fish of age (cm)

sex

sex code

m

transformed age for SFR parameterization of von Bertalanffy equation

n

sample size

Source

Kimura, D. K. 1990. Testing nonlinear regression paramters under heteroscedastic, normally distributed errors. Biometrics 46:697-708.

Description

Calculate the equilibrium Beverton-Holt estimator of instantaneous total mortality (Z) from length data with bootstrapped standard errors or the same using the Ehrhardt and Ault(1992) bias-correction

Usage

bheq(len, type = c(1,2), K = NULL, Linf = NULL, Lc = NULL, 
La = NULL, nboot = 100)

Arguments

Details

The standard Beverton-Holt equilibrium estimator of instantaneous total mortality (Z) from length data (page 365 in Quinn and Deriso (1999)) is calculated. The mean length for lengths >=Lc is calculated automatically. Missing data are removed prior to calculation. Estimates of standard error are made by bootstrapping length data >=Lc using package boot.

Value

Dataframe of length 1 containing mean length>=Lc, sample size>=Lc, Z estimate and standard error.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Ehrhardt, N. M. and J. S. Ault. 1992. Analysis of two length-based mortality models applied to bounded catch length frequencies. Trans. Am. Fish. Soc. 121:115-122.

Quinn, T. J. and R. B. Deriso. 1999. Quantitative Fish Dynamics. Oxford University Press, New York, New York. 542 pages.

See Also

bhnoneq

Examples

data(pinfish)
bheq(pinfish$sl,type=1,K=0.33,Linf=219.9,Lc=120,nboot=200)

Description

A nonequilibrium Beverton-Holt estimator of instantaneous total mortality (Z) from length data.

Usage

bhnoneq(year=NULL,mlen=NULL, ss=NULL, K = NULL, Linf = NULL, 
Lc = NULL, nbreaks = NULL, styrs = NULL, stZ = NULL, 
graph = TRUE)

Arguments

Details

The mean lengths for each year for lengths>=Lc. Following Gedamke and Hoening(2006), the model estimates nbreaks+1 Z values, the year(s) in which the changes in mortality began, the standard deviation of lengths>=Lc, and standard errors of all parameters. An AIC value is produced for model comparison. The estimated parameters for the number of nbreaks is equal to 2*nbreaks+2. Problematic parameter estimates may have extremely large t-values or extremely small standard error. Quang C. Huynh of Virginia Institute of Marine Science revised the function to make estimation more stable. Specifically, the derivative method BFGS is used in optim which allows more reliable convergence to the global minimum from a given set of starting values, a function is included to estimate Z assuming equilibrium, sigma is estimated analytically and convergence results . Use 0 nbreaks to get Z equilibrium.

Value

Note

Todd Gedamke provided the predicted mean length code in C++.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

Quang C. Huynh of Virginia Institute of Marine Science

References

Gedamke, T. and J. M. Hoenig. 2006. Estimating mortality from mean length data in nonequilibrium situations, with application to the assessment of goosefish. Trans. Am. Fish. Soc. 135:476-487

See Also

bheq

Examples

data(goosefish)
bhnoneq(year=goosefish$year,mlen=goosefish$mlen, ss=goosefish$ss,
K=0.108,Linf=126,Lc=30,nbreaks=1,styrs=c(1982),stZ=c(0.1,0.3))

Description

Growth increment data derived from tagging experiments on Pacific bonito (Sarda chiliensis) used to illustrate Francis's maximum likelihood method estimation of growth and growth variability (1988).

Usage

bonito

Format

A data frame with 138 observations on the following 4 variables.

T1

a numeric vector describing the release date

T2

a numeric vector describing the recovery date

L1

a numeric vector describing the length at release in cenitmeters

L2

a numeric vector describing the length at recapture in centimeters

Details

Note that Francis (1988) has discarded several records from the original dataset collected by Campbell et al. (1975).

Source

1

Francis, R.I.C.C., 1988. Maximum likelihood estimation of growth and growth variability from tagging data. New Zealand Journal of Marine and Freshwater Research, 22, p.42–51.

2

Campbell, G. & Collins, R., 1975. The age and growth of the Pacific bonito, Sarda chiliensis, in the eastern north Pacific. Calif. Dept. Fish Game, 61(4), p.181-200.

Description

Converts a log-mean and log-variance to the original scale and calculates confidence intervals

Usage

bt.log(meanlog = NULL, sdlog = NULL, n = NULL, alpha = 0.05)

Arguments

Details

There are two methods of calcuating the bias-corrected mean on the original scale. The bt.mean is calculated following equation 14 (the infinite series estimation) in Finney (1941). approx.bt.mean is calculated using the commonly known approximation from Finney (1941):

mean=exp(meanlog+sdlog^2/2). The variance is var=exp(2*meanlog)*(Gn(2*sdlog^2)-Gn((n-2)/(n-1)*sdlog^2) and standard deviation is sqrt(var) where Gn is the infinite series function (equation 10). The variance and standard deviation of the back-transformed mean are var.mean=var/n; sd.mean=sqrt(var.mean). The median is calculated as exp(meanlog). Confidence intervals for the back-transformed mean are from the Cox method (Zhou and Gao, 1997) modified by substituting the z distribution with the t distribution as recommended by Olsson (2005):

LCI=exp(meanlog+sdlog^2/2-t(df,1-alpha/2)*sqrt((sdlog^2/n)+(sdlog^4/(2*(n-1)))) and

UCI=exp(meanlog+sdlog^2/2+t(df,1-alpha/2)*sqrt((sdlog^2/n)+(sdlog^4/(2*(n-1))))

where df=n-1.

Value

A vector containing bt.mean, approx.bt.mean,var, sd, var.mean,sd.mean, median, LCI (lower confidence interval), and UCI (upper confidence interval).

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Finney, D. J. 1941. On the distribution of a variate whose logarithm is normally distributed. Journal of the Royal Statistical Society Supplement 7: 155-161.

Zhou, X-H., and Gao, S. 1997. Confidence intervals for the log-normal mean. Statistics in Medicine 16:783-790.

Olsson, F. 2005. Confidence intervals for the mean of a log-normal distribution. Journal of Statistics Education 13(1). www.amstat.org/publications/jse/v13n1/olsson.html

Examples

## The example below shows accuracy of the back-transformation
y<-rlnorm(100,meanlog=0.7,sdlog=0.2)
known<-unlist(list(known.mean=mean(y),var=var(y),sd=sd(y),
  var.mean=var(y)/length(y),sd.mean=sqrt(var(y)/length(y))))
est<-bt.log(meanlog=mean(log(y)),sdlog=sd(log(y)),n=length(y))[c(1,3,4,5,6)]
known;est

Description

The buffalo data frame has 20 rows and 3 columns. Cohort size and deaths for African buffalo from Sinclair (1977) as reported by Krebs (1989) in Table 12.1, page 415.

Usage

buffalo

Format

This data frame contains the following columns:

age

age interval

nx

number alive at start of each age interval

dx

number dying between age interval X and X+1

Source

Krebs, C. J. 1989. Ecological Methodologies. Harper and Row, New York, NY. 654 p.

Description

The catch data frame has 10 rows and 1 column.

Usage

catch

Format

This data frame contains the following columns:

value

catch data

Source

Massachusetts Division of Marine Fisheries

Description

Estimates selectivity-at-length from catch lengths and von Bertalanffy growth parameters.

Usage

catch.select(len = NULL, lenmin = NULL, binsize = NULL, 
peakplus = 1, Linf = NULL, K = NULL, t0 = NULL, subobs = FALSE)

Arguments

Details

This function applies the method of Pauly (1984) for calculating the selectivity-at-length from catch lengths and parameters from a von Bertalanffy growth curve. See Sparre and Venema(1998) for a detailed example of the application.

Length intervals are constructed based on the lenmin and binsize specified, and the maximum length observed in the data vector. Catch-at-length is tabularized using the lower and upper intervals and the data vector of lengths. The inclusion of a length in an interval is determined by lower interval>=length<upper interval. The age corresponding to the interval midpoint (t) is determined using the von Bertalanffy equation applied to the lower and upper bounds of each interval, summing the ages and dividing by 2. deltat is calculated for each interval using the equation: (1/k)*log((Linf-L1)/(Linf-L2)) where L1 and L2 are the lower and upper bounds of the length interval. log(catch/deltat) is the dependent variable and t is the predictor used in linear regression to estimate Z. Using the parameters from the catch curve analysis, "observed" selectivities (stobs) for the length intervals not included in the catch curve analysis are calculated using the equation: stobs=catch/(deltat*exp(a-Z*t)) where a and Z are the intercept and slope from the linear regression. The stobs values are transformed using an inverse logit (log(1/stobs-1)) and are regressed against t to obtain parameter estimates for the selectivity ogive. The estimated selectivity ogive (stest) is then calculated as stest=1/(1+exp(T1-T2*t)) where T1 and T2 are the intercept and slope from the log(1/stobs-1) versus t regression.

Value

list containing a dataframe with the lower and upper length intervals, the mid-point length interval, age corresponding to the interval mid-point, catch of the length interval, log(catch/deltat), the predicted log(catch/deltat) from the catch curve model fit (only for the peakplus interval and greater), the observed selectivities and the estimated selectivity, and two dataframes containing the parameters and their standard errors from the linear regressions for catch curve analysis and the selectivity ogive.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Pauly, D. 1984. Length-converted catch curves. A powerful tool for fisheries research in the tropics (Part III). ICLARM Fishbyte 2(1): 17-19.

Smith, M. W. and 5 others. 2012. Recommendations for catch-curve analysis. N. Am. J. Fish. Manage. 32: 956-967.

Sparre, P. and S. C. Venema. 1998. Introduction to tropical fish stock assessment. Part 1. Manual. FAO Fisheries Technical Paper, No. 206.1, Rev. 2. Rome. 407 p. Available on the world-wide web.

Examples

data(sblen)
catch.select(len=sblen$len_inches,binsize=2,lenmin=10,peakplus=1,Linf=47.5,K=0.15, 
t0=-0.3)

Description

This function estimates MSY following Martell and Froese(2012).

Usage

catchmsy(year = NULL, catch = NULL, catchCV = NULL, 
catargs = list(dist = "none", low = 0, up = Inf, unit = "MT"), 
l0 = list(low = 0, up = 1, step = 0), lt = list(low = 0, up = 1, 
refyr = NULL), 
sigv = 0, k = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), 
r = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), 
M = list(dist = "unif", low = 0.2, up = 0.2, mean = 0, sd = 0), 
nsims = 10000, catchout = 0, grout = 1, 
graphs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11), 
grargs = list(lwd = 1, pch = 16, cex = 1, nclasses = 20, mains = " ", 
cex.main = 1, 
cex.axis = 1, cex.lab = 1), 
pstats = list(ol = 1, mlty = 1, mlwd = 1.5, llty = 3, llwd = 1, ulty = 3, 
ulwd = 1), 
grtif = list(zoom = 4, width = 11, height = 13, pointsize = 10))

Arguments

Details

The method of Martell and Froese (2012) is used to produce estimates of MSY where only catch and information on resilience is known.

The Schaefer production model is

B[t+1]<-B[t]+r*B[t]*(1-B[t]/k)-catch[t]

where B is biomass in year t, r is the instrince rate of increase, k is the carrying capacity and catch is the catch in year t. Biomass is the first year is calculated by B[1]=k*l0. For sigv>0, the production equation is multiplied by exp(rnorm(1,0,sigv)) if process error is desired. The maximum sustainable yield (MSY) is calculated as

MSY=r*k/4

Biomass at MSY is calculated as

Bmsy=k/2

Fishing mortality at MSY is calculated as

Fmsy=r/2

Exploitation rate at MSY is calculated as

Umsy=(Fmsy/(Fmsy+M))*(1-exp(-Fmsy-M))

The overfishing limit in last year+1 is calculated as

OFL=B[last year +1]*Umsy

length(year)+1 biomass estimates are made for each run.

If using the R Gui (not Rstudio), run

graphics.off() windows(width=10, height=12,record=TRUE) .SavedPlots <- NULL

before running the catchmsy function to recall plots.

Value

The biomass estimates from each simulation are not stored in memory but are automatically saved to a .csv file named "Biotraj-cmsy.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values for each simulation (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims). This file is loaded to plot graph 11 and it must be present in the default or setwd() directory.

When catchout=1, catch values randomly selected are saved to a .csv file named "Catchtraj-cmsy.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims).

Use setwd() before running the function to change the directory where .csv files are stored.

Note

The random distribution function was adapted from Nadarajah, S and S. Kotz. 2006. R programs for computing truncated distributions. Journal of Statistical Software 16, code snippet 2.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Martell, S. and R. Froese. 2012. A simple method for estimating MSY from catch and resilience. Fish and Fisheries 14:504-514.

See Also

dbsra dlproj

Examples

## Not run: 
   data(lingcod)
   outpt<-catchmsy(year=lingcod$year,
     catch=lingcod$catch,catchCV=NULL,
     catargs=list(dist="none",low=0,up=Inf,unit="MT"),
    l0=list(low=0.8,up=0.8,step=0),
    lt=list(low=0.01,up=0.25,refyr=2002),sigv=0,
    k=list(dist="unif",low=4333,up=433300,mean=0,sd=0),
    r=list(dist="unif",low=0.015,up=0.1,mean=0,sd=0),
    M=list(dist="unif",low=0.18,up=0.18,mean=0.00,sd=0.00),
    nsims=30000)
 
## End(Not run)

Description

This function applies the catch-survey analysis method of Collie and Kruse (1998) for estimating abundance from catch and survey indices of relative abundance

Usage

catchsurvey(year = NULL, catch = NULL, recr = NULL, post = NULL, M = NULL,
 T = NULL, phi = NULL, w = 1, initial = c(NA,NA,NA),uprn = NA, graph = TRUE)

Arguments

Details

Details of the model are given in Collie and Kruse (1998).

Value

List containing the estimate of catchability, predicted recruit index by year (rest), estimate of recruit abundance (R), predicted post-recruit index by year (nest), post-recruit abundance (N), total abundance (TA: R+N), total instantaneous mortality (Z), and fishing mortality (Fmort)

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Collie JS and GH Kruse 1998. Estimating king crab (Paralithodes camtschaticus) abundance from commercial catch and research survey data. In: Jamieson GS, Campbell A, eds. Proceedings of the North Pacific Symposium on Invertebrate Stock Assessment and Management. Can Spec Publ Fish Aquat Sci. 125; p 73-83.

See also Collie JS and MP Sissenwine. 1983. Estimating population size from relative abundance data measured with error. Can J Fish Aquat Sci. 40:1871-1879.

Examples

## Example takes a bit of time to run
  ## Not run: 
   data(nshrimp)
   catchsurvey(year=nshrimp$year,catch=nshrimp$C,recr=nshrimp$r,post=nshrimp$n,M=0.25,
   T=0.5,phi=0.9,w=1,initial=c(500,500,0.7),uprn=10000)
## End(Not run)

Description

Statistical comparison of length frequencies is performed using the two-sample Kolmogorov & Smirnov test. Randomization procedures are used to derive the null probability distribution.

Usage

clus.lf(group = NULL, haul = NULL, len = NULL, number= NULL,
 binsize = NULL, resamples = 100)

Arguments

Details

Length frequency distributions of fishes are commonly tested for differences among groups (e.g., regions, sexes, etc.) using a two-sample Kolmogov-Smirnov test (K-S). Like most statistical tests, the K-S test requires that observations are collected at random and are independent of each other to satisfy assumptions. These basic assumptions are violated when gears (e.g., trawls, haul seines, gillnets, etc.) are used to sample fish because individuals are collected in clusters . In this case, the "haul", not the individual fish, is the primary sampling unit and statistical comparisons must take this into account.

To test for difference between length frequency distributions from simple random cluster sampling, a randomization test that uses "hauls" as the primary sampling unit can be used to generate the null probability distribution. In a randomization test, an observed test statistic is compared to an empirical probability density distribution of a test statistic under the null hypothesis of no difference. The observed test statistic used here is the Kolmogorov-Smirnov statistic (Ds) under a two-tailed test:

$Ds= max|S1(X)-S2(X)|$

where S1(X) and S2(X) are the observed cumulative length frequency distributions of group 1 and group 2 in the paired comparisons. S1(X) and S2(X) are calculated such that S(X)=K/n where K is the number of scores equal to or less than X and n is the total number of length observations (Seigel, 1956).

To generate the empirical probability density function (pdf), haul data are randomly assigned without replacement to the two groups with samples sizes equal to the original number of hauls in each group under comparison. The K-S statistic is calculated from the cumulative length frequency distributions of the two groups of randomized data. The randomization procedure is repeated resamples times to obtain the pdf of D. To estimate the significance of Ds, the proportion of all randomized D values that were greater than or equal to Ds is calculated.

It is assumed all fish caught are measured. If subsampling occurs, the number at length (measured) must be expanded to the total caught.

Data vectors described in arguments should be aggregated so that each record contains the number of fish in each length class by group and haul identifier. For example,

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Manly, B. F. J. 1997. Randomization, Bootstrap and Monte Carlos Methods in Biology. Chapman and Hall, New York, NY, 399 pp.

Seigel, S. 1956. Nonparametric Statistics for Behavioral Sciences. McGraw-Hill, New York, NY. 312 p.

See Also

clus.str.lf

Examples

data(codcluslen)
clus.lf(group=codcluslen$region,haul=c(paste("ST-",codcluslen$tow,sep="")),
 len=codcluslen$length, number=codcluslen$number, 
 binsize=5,resamples=100)

Description

Calculates mean attribute, variance, effective sample size, and degrees of freedom for samples collected by simple random cluster sampling.

Usage

clus.mean(popchar = NULL, cluster = NULL, clustotal = NULL, rho = NULL,
 nboot = 1000)

Arguments

Details

In fisheries, gears (e.g., trawls, haul seines, gillnets, etc.) are used to collect fishes. Often, estimates of mean population attributes (e.g., mean length) are desired. The samples of individual fish are not random samples, but cluster samples because the "haul" is the primary sampling unit. Correct estimation of mean attributes requires the use of cluster sampling formulae. Estimation of the general mean attribute and usual variance approximation follows Pennington et al. (2002). Variance of the mean is also estimated using the jackknife and bootstrap methods (Pennington and Volstad, 1994; Pennington et al., 2002). In addition, the effective sample size (the number of fish that would need to be sampled randomly to obtained the same precision as the mean estimate from cluster sampling) is also calculated for the three variance estimates. The total number of fish caught in a cluster (clustotal) allows correct computation for one- and two-stage sampling of individuals from each cluster (haul). In addition, if rho is specified, degrees of freedom are calculated by using Hedges (2007) for unequal cluster sizes (p. 166-167).

Value

Matrix table of total number of clusters (n), total number of samples (M), total number of samples measured (m), the mean attribute (R), usual variance approximation (varU), jackknife variance (varJ), bootstrap variance (varB), variance of population attribute (s2x), usual variance effective sample size (meffU), jackknife variance effective sample size, (meffJ), bootstrap variance effective sample size (meffB) and degrees of freedom (df) if applicable.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Hedges,L.V. 2007. Correcting a significance test for clustering. Journal of Educational and Behavioral Statistics. 32: 151-179.

Pennington, M. and J. J. Volstad. 1994. Assessing the effect of intra-haul correlation and variable density on estimates of population characteristics from marine surveys. Biometrics 50: 725-732.

Pennington, M. , L. Burmeister, and V. Hjellvik. 2002. Assessing the precision of frequency distributions estimated from trawl-survey samples. Fish. Bull. 100:74-80.

Examples

data(codcluslen)
temp<-codcluslen[codcluslen$region=="NorthCape" & codcluslen$number>0,]
temp$station<-c(paste(temp$region,"-",temp$tow,sep=""))
total<-aggregate(temp$number,list(temp$station),sum)
names(total)<-c("station","total")
temp<-merge(temp,total,by.x="station",by.y="station")
newdata<-data.frame(NULL)
for(i in 1:as.numeric(length(temp[,1]))){
  for(j in 1:temp$number[i]){
    newdata<-rbind(newdata,temp[i,])
  }
}
clus.mean(popchar=newdata$length,cluster=newdata$station,
         clustotal=newdata$total)

Description

Calculates the intracluster correlation coefficients according to Lohr (1999) and Donner (1986) for a single group

Usage

clus.rho(popchar=NULL, cluster = NULL, type = c(1,2,3), est = 0, nboot = 500)

Arguments

Details

The intracluster correlation coefficient (rho) provides a measure of similarity within clusters. type = 1 is defined to be the Pearson correlation coefficient for NM(M-1)pairs (yij,yik) for i between 1 and N and j<>k (see Lohr (1999: p. 139). The average cluster size is used as the equal cluster size quantity in Equation 5.8 of Lohr (1999). If the clusters are perfectly homogeneous (total variation is all between-cluster variability), then ICC=1.

type = 2 is the adjusted r-square, an alternative quantity following Equation 5.10 in Lohr (1999). It is the relative amount of variability in the population explained by the cluster means, adjusted for the number of degrees of freedom. If the clusters are homogeneous, then the cluster means are highly variable relative to variation within clusters, and the r-square will be high.

type = 3 is calculated using one-way random effects models (Donner, 1986). The formula is

rho = (BMS-WMS)/(BMS+(m-1)*WMS)

where BMS is the mean square between clusters, WMS is the mean square within clusters and m is the adjusted mean cluster size for clusters with unequal sample size. All clusters with zero elementary units should be deleted before calculation. type = 3 can be used with binary data (Ridout et al. 1999)

If est=1, the boostrap mean (value), variance of the mean and 0.025 and 0.975 percentiles are outputted.

Value

rho values, associated statistics, and bootstrap replicates

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Donner, A. 1986. A review of inference procedures for the intraclass correlation coefficient in the one-way random effects model. International Statistical Review. 54: 67-82.

Lohr, S. L. Sampling: design and analysis. Duxbury Press,New York, NY. 494 p.

Ridout, M. S., C. G. B. Demetrio, and D. Firth. 1999. Estimating intraclass correlation for binary data. Biometrics 55: 137-148.

See Also

clus.lf clus.str.lf clus.mean

Examples

data(codcluslen)
  tem<-codcluslen[codcluslen[,1]=="NorthCape" & codcluslen[,3]>0,]
  outs<-data.frame(tow=NA,len=NA)
  cnt<-0
  for(i in 1:as.numeric(length(tem$number))){
    for(j in 1:tem$number[i]){
     cnt<-cnt+1
     outs[cnt,1]<-tem$tow[i]
     outs[cnt,2]<-tem$length[i]
   }
 }
 clus.rho(popchar=outs$len,cluster=outs$tow)

Description

Calculates a common intracluster correlation coefficients according to Donner (1986: 77-79) for two or more groups with unequal cluster sizes, and tests for homogeneity of residual error among groups and a common coefficient among groups.

Usage

clus.rho.g(popchar=NULL, cluster = NULL, group = NULL)

Arguments

Details

The intracluster correlation coefficient (rho) provides a measure of similarity within clusters. rho is calculated using a one-way nested random effects model (Donner, 1986: 77-79). The formula is

rho = (BMS-WMS)/(BMS+(m-1)*WMS)

where BMS is the mean square among clusters within groups, WMS is the mean square within clusters and m is the adjusted mean cluster size for clusters with unequal sample sizes. All clusters with zero elementary units should be deleted before calculation. In addition, approximate 95 are generated and a significance test is performed.

Bartlett's test is used to determine if mean square errors are constant among groups. If Bartlett's test is not significant, the test for a common correlation coefficient among groups is valid.

Value

rho value and associate statistics

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Donner, A. 1986. A review of inference procedures for the intraclass correlation coefficient in the one-way random effects model. International Statistical Review. 54: 67-82.

See Also

clus.str.lf clus.lf clus.mean

Examples

data(codcluslen)
   temp<-codcluslen[codcluslen$number>0,]
   temp$station<-c(paste(temp$region,"-",temp$tow,sep=""))
   total<-aggregate(temp$number,list(temp$station),sum)
   names(total)<-c("station","total")
   temp<-merge(temp,total,by.x="station",by.y="station")
   newdata<-data.frame(NULL)
   for(i in 1:as.numeric(length(temp[,1]))){
    for(j in 1:temp$number[i]){
     newdata<-rbind(newdata,temp[i,])
    }
  }
  newdata<-newdata[,-c(5)]
 clus.rho.g(popchar=newdata$length,cluster=newdata$station,group=newdata$region)

Description

Statistical comparison of length frequencies is performed using the two-sample Kolmogorov & Smirnov test. Randomization procedures are used to derive the null probability distribution.

Usage

clus.str.lf(group = NULL, strata = NULL, weights = NULL,
 haul = NULL, len = NULL, number = NULL, binsize = NULL,
 resamples = 100)

Arguments

Details

Length frequency distributions of fishes are commonly tested for differences among groups (e.g., regions, sexes, etc.) using a two-sample Kolmogov-Smirnov test (K-S). Like most statistical tests, the K-S test requires that observations are collected at random and are independent of each other to satisfy assumptions. These basic assumptions are violated when gears (e.g., trawls, haul seines, gillnets, etc.) are used to sample fish because individuals are collected in clusters . In this case, the "haul", not the individual fish, is the primary sampling unit and statistical comparisons must take this into account.

To test for difference between length frequency distributions from stratified random cluster sampling, a randomization test that uses "hauls" as the primary sampling unit can be used to generate the null probability distribution. In a randomization test, an observed test statistic is compared to an empirical probability density distribution of a test statistic under the null hypothesis of no difference. The observed test statistic used here is the Kolmogorov-Smirnov statistic (Ds) under a two-tailed test:

$Ds= max|S1(X)-S2(X)|$

where S1(X) and S2(X) are the observed cumulative proportions at length for group 1 and group 2 in the paired comparisons.

Proportion of fish of length class j in strata-set (group variable) used to derive Ds is calculated as

$p_j=\frac{\sum{A_k\bar{X}}_{jk}}{\sum{A_k\bar{X}}_k}$

where $A_k$ is the weight of stratum k, $\bar{X}_{jk}$ is the mean number per haul of length class j in stratum k, and $\bar{X}_k$ is the mean number per haul in stratum k. The numerator and denominator are summed over all k. Before calculation of cumulative proportions, the length class distributions for each group are corrected for missing lengths and are constructed so that the range and intervals of each distribution match.

It is assumed all fish caught are measured. If subsampling occurs, the numbers at length (measured) must be expanded to the total caught.

To generate the empirical probability density function (pdf), length data of hauls from all strata are pooled and then hauls are randomly assigned without replacement to each stratum with haul sizes equal to the original number of stratum hauls. Cumulative proportions are then calculated as described above. The K-S statistic is calculated from the cumulative length frequency distributions of the two groups of randomized data. The randomization procedure is repeated resamples times to obtain the pdf of D. To estimate the significance of Ds, the proportion of all randomized D values that were greater than or equal to Ds is calculated (Manly, 1997).

Data vectors described in arguments should be aggregated so that each record contains the number of fish in each length class by group, strata, weights, and haul identifier. For example,

To correctly calculate the stratified mean number per haul, zero tows must be included in the dataset. To designate records for zero tows, fill the length class and number at length with zeros. The first line in the following table shows the appropriate coding for zero tows:

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Manly, B. F. J. 1997. Randomization, Bootstrap and Monte Carlos Methods in Biology. Chapman and Hall, New York, NY, 399 pp.

Seigel, S. 1956. Nonparametric Statistics for Behavioral Sciences. McGraw-Hill, New York, NY. 312 p.

See Also

clus.lf

Examples

data(codstrcluslen)
clus.str.lf(
group=codstrcluslen$region,strata=codstrcluslen$stratum,
 weights=codstrcluslen$weights,haul=codstrcluslen$tow,
 len=codstrcluslen$length,number=codstrcluslen$number,
 binsize=5,resamples=100)

Description

Calculates Hedges (2007) t-statistic adjustment and degrees of freedom for a t-test assuming unequal variances and clustered data with clusters of unequal size.

Usage

clus.t.test(popchar = NULL, cluster = NULL, group = NULL,
      rho = NULL, alpha = 0.05, alternative = c("two.sided"))

Arguments

Details

A two-sample t-test with unequal variances (Sokal and Rohlf, 1995) is performed on clustered data. The t-statistic and degrees of freedom are corrected for clustering according to Hedges (2007).

Value

List with null hypothesis of test and matrix table with mean of each group, rho, ntilda (Equation 14 of Hedges 2007), nu (Equation 15), degrees of freedom (Equation 16), uncorrected t-statistic, cu (Equation 18), the t-statistic adjusted for clustering, critical t value for degrees of freedom and alpha, and probability of significance.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Sokal,R.R.and F.J.Rohlf. 1995. Biometry, 3rd Edition. W.H. Freeman and Company, New York, NY. 887 p.

Hedges,L.V. 2007. Correcting a significance test for clustering. Journal of Educational and Behavioral Statistics. 32: 151-179.

Examples

data(codcluslen)
   temp<-codcluslen[codcluslen$number>0,]
   temp$station<-c(paste(temp$region,"-",temp$tow,sep=""))
   total<-aggregate(temp$number,list(temp$station),sum)
   names(total)<-c("station","total")
   temp<-merge(temp,total,by.x="station",by.y="station")
   newdata<-data.frame(NULL)
   for(i in 1:as.numeric(length(temp[,1]))){
    for(j in 1:temp$number[i]){
     newdata<-rbind(newdata,temp[i,])
    }
  }
 newdata<-newdata[,-c(5)]
 clus.t.test(popchar=newdata$length,cluster=newdata$station,
            group=newdata$region,rho=0.72,
            alpha=0.05,alternative="two.sided")

Description

Fits the von Bertalanffy growth equation to clustered length and age by using nonlinear least-squares and by bootstrapping clusters

Usage

clus.vb.fit(len = NULL, age = NULL, cluster = NULL, nboot = 1000,
sumtype = 1, control = list(maxiter=10000, minFactor=1/1024,tol=1e-5))

Arguments

Details

A standard von Bertalanffy growth curve is fitted to length-at-age data of each nboot sample of clusters using nonlinear least-squares (function nls). Standard errors are calculated using function sd applied to bootstrap parameters.

Value

List containing a summary of successful model fits and parameter estimates, standard errors and 95 percent confidence intervals, and the average correlation matrix.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

Examples

## Not run: 
	data(pinfish)
	with(pinfish,clus.vb.fit(len=sl,age=age,cluster=field_no,nboot=100))
  
## End(Not run)

Description

The codcluslen data frame has 334 rows and 4 columns.

Usage

codcluslen

Format

This data frame contains the following columns:

region

NorthCape = North of Cape Cod; SouthCape =South of Cape Cod

tow

Tow number

length

Length class (total length, cm)

number

Number in length class

Source

Massachusetts Division of Marine Fisheries

Description

The codstrcluslen data frame has 334 rows and 6 columns.

Usage

codstrcluslen

Format

This data frame contains the following columns:

region

NorthCape = North of Cape Cod; SouthCape = South of Cape Cod

stratum

Stratum number

tow

Tow number

weights

Stratum area (square nautical-miles)

length

Length class (total length cm)

number

Number in length class

Source

Massachusetts Division of Marine Fisheries, 30 Emerson Avenue, Gloucester, MA 01930

Description

This function takes multiple mean and sample variance estimates and combines them.

Usage

combinevar(xbar = NULL, s_squared = NULL, n = NULL)

Arguments

Details

If a Monte Carlo simulation is run over 1000 loops and then again over another 1000 loops, one may wish to update the mean and variance from the first 1000 loops with the second set of simulation results.

Value

Vector containing the combined mean and sample variance.

Author(s)

John M. Hoenig, Virginia Institute of Marine Science [email protected]

Examples

xbar <- c(5,5)
s<-c(2,4)
n <- c(10,10)
combinevar(xbar,s,n)

Description

Compute likelihood ratio tests for two or more growthlrt.plus model objects via Kimura (1990)

Usage

compare.lrt.plus(...)

Arguments

Details

Likelihood ratio and F tests are computed for models compared against one another in the order specified.

Value

List containing model test statistics

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Kimura, D. K. 1990. Testing nonlinear reression parameters under heteroscedastic, normally distributed errors. Biometrics 46: 697-708.

See Also

growthlrt.plus

Examples

## This is a typical specification, not a working example
## Not run: 
compare.lrt.plus(model1,model2)
## End(Not run)

Description

Function compares graphically the readings of two age readers and calculates 2 chi-square statistics for tests of symmetry.

Usage

compare2(readings, usecols = c(1,2), twovsone = TRUE, plot.summary = TRUE, 
  barplot = TRUE, chi = TRUE, pool.criterion = 1, cont.cor = TRUE, 
  correct = "Yates", first.name = "Reader A",second.name = "Reader B")

Arguments

Details

This function can plot the number of readings of age j by reader 2 versus the number of readings of age i by reader 1 (if twovsone=TRUE). Optionally, it will add the number of readings above, on, and below the 45 degree line (if plot.summary=TRUE). The function can make a histogram of the differences in readings (if barplot=TRUE). Finally, the program can calculate 2 chi-square test statistics for tests of the null hypothesis that the two readers are interchangeable vs the alternative that there are systematic differences between readers (if chi=TRUE). The tests are tests of symmetry (Evans and Hoenig, 1998). If cont.cor=T, then correction for continuity is applied to the McNemar-like chi-square test statistic; the default is to apply the Yates correction but if correct="Edwards" is specified then the correction for continuity is 1.0 instead of 0.5.

Value

Separate lists with tables of various statistics associated with the method.

Author(s)

John Hoenig, Virginia Institute of Marine Science, 18 December 2012. [email protected]

References

Evans, G.T. and J.M. Hoenig. 1998. Viewing and Testing Symmetry in Contingency Tables, with Application to Fish Ages. Biometrics 54:620-629.).

Examples

data(sbotos)
 compare2(readings=sbotos,usecols=c(1,2),twovsone=TRUE,plot.summary=TRUE,
 barplot=FALSE,chi=TRUE,pool.criterion=1,cont.cor=TRUE,correct="Yates",
 first.name="Reader A",second.name="Reader B")

Description

Convert instantaneous fishing mortality rate (F) to annual exploitation rate (mu) and vice versa for Type I and II fisheries.

Usage

convmort(value = NULL, fromto = 1, type = 2, M = NULL)

Arguments

Details

Equations 1.6 and 1.11 of Ricker (1975) are used.

Value

A vector of the same length as value containing the converted values.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Ricker, W. E. 1975. Computation and interpretation of biological statistics of fish populations. Bull. Fish. Res. Board. Can. 191: 382 p.

Examples

convmort(0.3,fromto=1,type=2,M=0.15)

Description

The counts data frame has 31 rows and 2 columns. Run size data of alewife (Alosa pseudoharengus) in Herring River, MA from 1980-2010

Usage

counts

Format

This data frame contains the following columns:

year

vector of run year

number

vector of run counts in number of fish

Source

Massachusetts Division of Marine Fisheries, 30 Emerson Avenue, Gloucester, MA

Description

Cowcod catch data from literature sources in Martell and Froese (2012).

Usage

cowcod

Format

A data frame with 109 observations on the following 2 variables.

year

a numeric vector describing the year of catch

catch

a numeric vector describing the annual catch in metric tons

Description

Calculates the mean cpue after replacing unusually large catches with expected values using the method of Kappenman (1999)

Usage

cpuekapp(x = NULL, nlarge = NULL, absdif = 0.001)

Arguments

Details

Use function gap to choose the number of unusually large values.

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Kappenman, R. F. 1999. Trawl survey based abundance estimation using data sets with unusually large catches. ICES Journal of Marine Science. 56: 28-35.

See Also

gap

Examples

## Not run: 
  ## Data from Table 1 in Kappenman (1999)
  data(kappenman)
  cpuekapp(kappenman$cpue,1)
 
## End(Not run)

Description

The darter data frame has 7 rows and 2 columns. Sequence of catch data for the faintail darter from removal experiments by Mahon as reported by White et al.(1982). This dataset is often use to test new depletion estimators because the actual abundance is known (N=1151).

Usage

darter

Format

This data frame contains the following columns:

catch

catch data

effort

constant effort data

Source

White, G. C., D. R. Anderson, K. P. Burnham, and D. L. Otis. 1982. Capture-recapture and Removal Methods for Sampling Closed Populations. Los Alamos National Laboratory LA-8787-NERP. 235 p.

Description

This function estimates MSY from catch following Dick and MAcCall (2011).

Usage

dbsra(year = NULL, catch = NULL, catchCV = NULL, 
catargs = list(dist = "none", low = 0, up = Inf, unit = "MT"), 
agemat = NULL, maxn=25, k = list(low = 0, up = NULL, tol = 0.01, permax = 1000), 
b1k = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0),
btk = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0, refyr = NULL), 
fmsym = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), 
bmsyk = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), 
M = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), nsims = 10000, 
catchout = 0, grout = 1, 
graphs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15), 
grargs = list(lwd = 1, cex = 1, nclasses = 20, mains = " ", cex.main = 1, 
cex.axis = 1, 
cex.lab = 1), pstats = list(ol = 1, mlty = 1, mlwd = 1.5, llty = 3, llwd = 1, 
ulty = 3, ulwd = 1), 
grtif = list(zoom = 4, width = 11, height = 13, pointsize = 10))

Arguments

Details

The method of Dick and MAcCall (2011) is used to produce estimates of MSY where only catch and information on resilience and current relative depletion is known.

The delay-difference model is used to propogate biomass:

B[t+1]<-B[t]+P[Bt-a]-C[t]

where B[t] is biomass in year t, P[Bt-a] is latent annual production based on parental biomass agemat years earlier and C[t] is the catch in year t. Biomass in the first year is assumed equal to k.

If Bmsy/k>=0.5, then P[t] is calculated as

P[t]<-g*MSY*(B[t-agemat]/k)-g*MSY*(B[t-agemat]/k)^n

where MSY is k*Bmsy/k*Umsy, n is solved iteratively using the equation, Bmsy/k=n^(1/(1-n)), and g is (n^(n/(n-1)))/(n-1). Fmsy is calculated as Fmsy=Fmsy/M*M and Umsy is calculated as (Fmsy/(Fmsy+M))*(1-exp(-Fmsy-M)).

If Bsmy/k < 0.5, Bjoin is calculated based on linear rules: If Bmsy/k<0.3, Bjoin=0.5*Bmsy/k*k If Bmsy/k>0.3 and Bmsy/k<0.5, Bjoin=(0.75*Bmsy/k-0.075)*k

If any B[t-a]<Bjoin, then the Schaefer model is used to calculated P:

P[Bt-agematt<Bjoin]<-B[t-agemat]*(P(Bjoin)/Bjoin+c(B[t-agemat]-Bjoin))

where c =(1-n)*g*MSY*Bjoin^(n-2)*K^(-n)

Biomass at MSY is calculated as: Bmsy=(Bmsy/k)*k

The overfishing limit (OFL) is Umsy*B[termyr].

length(year)+1 biomass estimates are made for each run.

The rule for accepting a run is: if(min(B)>0 && max(B)<=k &&

(objective function minimum<=tol^2) && abs(((max(B)-k)/k)*100)<=permax && n<=maxn

If using the R Gui (not Rstudio), run

graphics.off() windows(width=10, height=12,record=TRUE) .SavedPlots <- NULL

before running the dbsra function to recall plots.

Value

The biomass estimates from each simulation are not stored in memory but are automatically saved to a .csv file named "Biotraj-dbsra.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values for each simulation (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims). This file is loaded to plot graph 13 and it must be present in the default or setwd() directory.

When catchout=1, catch values randomly selected are saved to a .csv file named "Catchtraj-dbsra.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims).

Use setwd() before running the function to change the directory where .csv files are stored.

Note

The random distribution function was adapted from Nadarajah, S and S. Kotz. 2006. R programs for computing truncated distributions. Journal of Statistical Software 16, code snippet 2.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Dick, E. J. and A. D. MacCall. 2011. Depletion-based stock reduction analysis: a catch-based method for determining sustainable yield for data-poor fish stocks. Fisheries Research 110: 331-341.

See Also

catchmsy dlproj

Examples

## Not run: 
  data(cowcod)
  dbsra(year =cowcod$year, catch = cowcod$catch, catchCV = NULL, 
    catargs = list(dist="none",low=0,up=Inf,unit="MT"),
    agemat=11, k = list(low=100,up=15000,tol=0.01,permax=1000), 
    b1k = list(dist="none",low=0.01,up=0.99,mean=1,sd=0.1),
    btk = list(dist="beta",low=0.01,up=0.99,mean=0.4,sd=0.1,refyr=2009),
    fmsym = list(dist="lnorm",low=0.1,up=2,mean=-0.223,sd=0.2),
    bmsyk = list(dist="beta",low=0.05,up=0.95,mean=0.4,sd=0.05),
    M = list(dist="lnorm",low=0.001,up=1,mean=-2.90,sd=0.4),
    nsims = 10000)
 
## End(Not run)

Description

Calculates the mean and variance of a catch series based on the delta distribution described in Pennington (1983).

Usage

deltadist(x = NULL)

Arguments

Details

Data from marine resources surveys usually contain a large proportion of hauls with no catches. Use of the delta-distribution can lead to more efficient estimators of the mean and variance because zeros are treated separately. The methods used here to calculate the delta distribution mean and variance are given in Pennington (1983).

Value

vector containing the delta mean and associated variance.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Pennington, M. 1983. Efficient estimators of abundance for fish and plankton surveys. Biometrics 39: 281-286.

Examples

data(catch)
deltadist(catch$value)

Description

Variable and constant effort models for the estimation of abundance from catch-effort depletion data assuming a closed population.

Usage

deplet(catch = NULL, effort = NULL, method = c("l", "d", "ml",
 "hosc", "hesc", "hemqle", "wh"), kwh=NULL, nboot = 500, Nstart=NULL)

Arguments

Details

The variable effort models include the Leslie-Davis (l) estimator (Leslie and Davis, 1939), the effort-corrected Delury (d) estimator (Delury,1947; Braaten, 1969), the maximum likelihood (ml) method of Gould and Pollock (1997), sample coverage estimator for the homogeneous model (hosc) of Chao and Chang (1999), sample coverage estimator for the heterogeneous model (hesc) of Chao and Chang (1999), and the maximum quasi-likelihood estimator for the heterogeneous model (hemqle) of Chao and Chang (1999). The variable effort models can be applied to constant effort data by simply filling the effort vector with 1s. Three removals are required to use the Leslie, Delury, and Gould and Pollock methods.

The constant effort model is the generalized removal method of Otis et al. 1978 reviewed in White et al. (1982: 109-114). If only two removals, the two-pass estimator of N in White et al. (1982:105) and the variance estimator of Otis et al. (1978: 108) are used.

Note: Calculation of the standard error using the ml method may take considerable time.

For the Delury method, zero catch values are not allowed because the log-transform is used.

For the generalized removal models, if standard errors appear as NAs but parameter estimates are provided, the inversion of the Hessian failed. If parameter estimates and standard errors appear as NAs, then model fitting failed.

For the Chao and Chang models, if the last catch value is zero, it is deleted from the data. Zero values between positive values are permitted.

Value

Separate output lists with the method name and extension .out are created for each method and contain tables of various statistics associated with the method.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Braaten, D. O. 1969. Robustness of the Delury population estimator. J. Fish. Res. Board Can. 26: 339-355.

Chao, A. and S. Chang. 1999. An estimating function approach to the inference of catch-effort models. Environ. Ecol. Stat. 6: 313-334.

Delury, D. B. 1947. On the estimation of biological populations. Biometrics 3: 145-167.

Gould, W. R. and K. H. Pollock. 1997. Catch-effort maximum likelihood estimation of important population parameters. Can. J. Fish. Aquat. Sci 54: 890-897.

Leslie, P. H. and D. H.S. Davis. 1939. An attempt to determine the absolute number of rats on a given area. J. Anim. Ecol. 9: 94-113.

Otis, D. L., K. P. Burnham, G. C. White, and D. R. Anderson. 1978. Statistical inference from capture data on closed animal populations. Wildl. Monogr. 62: 1-135.

White, G. C., D. R. Anderson, K. P. Burnham, and D. L. Otis. 1982. Capture-recapture and Removal Methods for Sampling Closed Populations. Los Alamos National Laboratory LA-8787-NERP. 235 p.

Examples

data(darter)
deplet(catch=darter$catch,effort=darter$effort,method="hosc") 
hosc.out

Description

Make biomass projections by using inputted catch and results of dbsra or catchmsy functions

Usage

dlproj(dlobj = NULL, projyears = NULL, projtype = 1, projcatch = NULL, 
grout = 1, grargs = list(lwd = 1, unit = "MT", mains = " ", cex.main = 1, 
cex.axis = 1, cex.lab = 1), grtif = list(zoom = 4, width = 11, height = 13, 
pointsize = 10))

Arguments

Details

The biomass estimate of the last year+1 is used as the starting biomass (year 1 in projections) and leading parameters from each plausible (accepted) run are used to project biomass ahead projyears years using either the MSY estimate (median or mean) from all plausible runs or inputted catch values. The biomass estimates are loaded from either the "Biotraj-dbsra.csv" or "Biotroj-cmsy.csv" files that were automatically saved in functions "dbsra" and "catchmsy".

Use setwd() before running the function to change the directory where .csv files are stored.

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Martell, S. and R. Froese. 2012. A simple method for estimating MSY from catch and resilience. Fish and Fisheries 14:504-514.

Dick, E. J. and A. D. MacCall. 2011. Depletion-based stock reduction analysis: a catch-based method for determining sustainable yield for data-poor fish stocks. Fisheries Research 110: 331-341.

See Also

catchmsy dbsra

Examples

## Not run: 
  data(lingcod)
   outs<-catchmsy(year=lingcod$year,
    catch=lingcod$catch,catchCV=NULL,
    catargs=list(dist="none",low=0,up=Inf,unit="MT"),
    l0=list(low=0.8,up=0.8,step=0),
    lt=list(low=0.01,up=0.25,refyr=2002),sigv=0,
    k=list(dist="unif",low=4333,up=433300,mean=0,sd=0),
    r=list(dist="unif",low=0.015,up=0.5,mean=0,sd=0),
    bk=list(dist="unif",low=0.5,up=0.5,mean=0,sd=0),
    M=list(dist="unif",low=0.24,up=0.24,mean=0.00,sd=0.00),
    nsims=30000)
   outbio<-dlproj(dlobj = outs, projyears = 20, projtype = 0, grout = 1)
 
## End(Not run)

Description

Estimation of von Bertanffy growth parameters based on length-stratified age samples (Perrault et al., 2020)

Usage

ep_growth(len=NULL,age=NULL,Nh=NULL,nh=NULL,starts=list(Linf=60,
k=0.1,a0=-0.01,CV=0.5), 
bin_size=2,nlminb.control=list(eval.max=5000, 
iter.max=5000,trace=10),
tmb.control=list(maxit=5000,trace=FALSE),plot=TRUE)

Arguments

Details

The von Bertalanffy growth model Lage=Linf*(1-exp(-K*(age-a0)) is fitted to length-at-age data collected via length-stratified sampling following the EP method of Perreault et al. (2020). A plot of model fit is generated unless plot=FALSE.

Value

List containing list elements of the model convergence, parameter estimates and predicted values.

Author(s)

Andrea Perrault, Marine Institute of Memorial University of Newfoundland

[email protected]

References

Perrault, A. M. J., N. Zhang and Noel G. Cadigan. 2020. Estimation of growth parameters based on length-stratified age samples. Canadian Journal of Fisheries and Aquatic Sciences 77: 439-450.

Examples

## Not run: 
	data(ep.data)
	ep_growth(len=ep.data$length, age=ep.data$age, Nh=ep.data$Nh, nh=ep.data$nh,
        starts=list(Linf=60,
        k=0.1, CV=0.5 ,a0=-0.01), bin_size=2,
          nlminb.control=list(eval.max=5000, iter.max=5000, trace=10), 
          tmb.control=list(maxit=5000, trace=F),plot=TRUE)
 
## End(Not run)

Description

The catch data frame has 1072 rows and 4 columns.

Usage

ep.data

Format

This data frame contains the following columns:

length

length in cm.

age

age of fish (yrs).

Nh

the total sample size per bin (includes unaged fish).

nh

the total aged sample size per bin.

Source

Andrea Perrault, Marine Institute of Memorial University of Newfoundland

Description

Eggs-per-recruit(EPR) analysis is conducted following Gabriel et al. (1989) except fecundity-at-age is substituted for weight-at-age. Reference points of F and EPR for percentage of maximum spawning potential are calculated.

Usage

epr(age = NULL, fecund = NULL, partial = NULL, pmat = pmat,
 M = NULL, pF = NULL, pM = NULL, MSP = 40, plus = FALSE,
 oldest = NULL, maxF = 2, incrF = 1e-04)

Arguments

Details

Eggs-per-recruit analysis is conducted following Gabriel et al. (1989). The F and EPR for the percentage maximum spawning potential reference point are calculated. If the last age is a plus-group, the cohort is expanded to the oldest age and the fecund, partial, pmat, and M values for the plus age are applied to the expanded cohort ages.

Value

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Gabriel, W. L., M. P. Sissenwine, and W. J. Overholtz. 1989. Analysis of spawning stock biomass per recruit: an example for Georges Bank haddock. North American Journal of Fisheries Management 9: 383-391.

See Also

ypr sbpr

Examples

data(menhaden)
	epr(age=menhaden$age,fecund=menhaden$fecundity,partial=menhaden$partial,
	pmat=menhaden$pmat,M=menhaden$M,pF=0,pM=0,MSP=40,plus=TRUE,maxF=4,incrF=0.01,oldest=10)

Description

Check design of parameter structure before use in function fm_telemetry.

Usage

fm_checkdesign(occasions = NULL, design = NULL, type = "F" )

Arguments

Details

The program allows the configuration of different parameter structure for the estimation of fishing and natural mortalities, and detection probabilities. These structures are specified in design. Consider the following examples:

Example 1

Tags are relocated over seven occasions. One model structure might be constant fishing mortality estimates over occasions 1-3 and 4-6. To specify this model structure: design is c(“1”,“4”).

Note: The structures of design must always contain the first occasion for fishing mortality and natural mortality, whereas the structure for the probability of detection must not contain the first occasion.

Example 2

Tags are relocated over six occasions. One model structure might be separate fishing mortality estimates for occasion 1-3 and the same parameter estimates for occasions 4-6. The design is c(“1:3*4:6”).

Note: The structures of Fdesign and Mdesign must always start with the first occasion, whereas the structure for Pdesign must always start with the second occasion.

Use the multiplication sign to specify occasions whose estimates of F, M or P will be taken from values of other occasions.

Example 3

Specification of model 3 listed in Table 1 of Hightower et al. (2001) is shown. Each occasion represented a quarter of the year. The quarter design for F specifies that quarterly estimates are the same in both years. design is c(“1*14”,“4*17”,“7*20”,“11*24”).

Example 4

In Hightower et al. (2001), the quarter and year design specifies that estimates are made for each quarter but are different for each year. design is

c(“1”, “4”, “7”, “11”, “14”, “17”, “20”, “24”).

If the number of occasions to be assigned parameters from other occasions are less than the number of original parameters (e.g., c(“11:13*24:25”), then only the beginning sequence of original parameters equal to the number of occasions are used. For instance, in c(“11:13*24:25”), only parameters 11 and 12 would be assigned to occasions 24 and 25.

If the number of occasions to be assigned parameters from other occasions are greater than the number of original parameters (e.g., c(“11:12*24:26”)), then the last original parameter is re-cycled. In the example c(“11:12*24:26”), the parameter for occasion 12 is assigned to occasions 25 and 26.

Value

dataframe containing the parameter order by occasion.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

See Also

fm_telemetry

Examples

fm_checkdesign(occasions=27, design=c("1*14","4*17","7*20","11*24"),type="F")

Description

Calculates model averaged estimates of instantaneous fishing, natural and probability of detection for telemetry models of Hightower et al. (2001).

Usage

fm_model_avg(..., global = NULL, chat = 1)

Arguments

Details

Model estimates are generated from function fm_telemetry. Averaging of model estimates follows the procedures in Burnham and Anderson (2002). Variances of parameters are adjusted for overdispersion using the c-hat estimate from the global model : sqrt(var*c-hat). If c-hat of the global model is <1, then c-hat is set to 1. The c-hat is used to calculate the quasi-likelihood AIC and AICc metrics for each model (see page 69 in Burnham and Anderson(2002)). QAICc differences among models are calculated by subtracting the QAICc of each model from the model with the smallest QAICc value. These differences are used to calculate the Akaike weights for each model following the formula on page 75 of Burnham and Anderson (2002). The Akaike weights are used to calculate the weighted average and standard error of parameter estimates by summing the product of the model-specific Akaike weight and parameter estimate across all models. An unconditional standard error is also calculated by sqrt(sum(QAICc wgt of model i * (var of est of model i + (est of model i - avg of all est)^2))).

Value

List containing model summary statistics, model-averaged estimates of fishing, natural and probability of detections and their weighted and uncondtional standard errors .

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Burnham, K. P. and D. R. Anderson. 2002. Model selection and multimodel inference : A Practical Information-Theorectic Approach, 2nd edition. Spriner-Verlag, New York, NY. 488 p.

See Also

fm_telemetry

Examples

## This is a typical specification, not a working example
## Not run: 
fm_model_avg(model1,model2,model3,model4,model5,model6,model7,global="model7")
## End(Not run)

Description

The method of Hightower et al. (2001) is implemented to estimate fishing mortality, natural mortality and probability of detection from telemetry data.

Usage

fm_telemetry(filetype = c(1), caphistory = NULL, Fdesign = NULL, Mdesign = NULL, 
Pdesign = NULL, whichlivecells =  NULL, 
whichdeadcells = NULL, constant = 1e-14, initial = NULL, 
invtol = 1e-44, control = list(reltol=1e-8,maxit=1000000))

Arguments

Details

The telemetry method of Hightower et al. (2001) is implemented. Individual capture histories are used in the function. The function uses complete capture histories (Burnham et al., 1987) and it is the presence of specific codes in the individual capture histories that split the capture histories into live and dead arrays. F and M estimates are needed for occasions 1 to the total number of occasions minus 1 and P estimates are needed for occasions 2 to the total number of occasions.

Capture histories are coded following Burnham et al. (1987)(i.e., 0 = not relocated, and 1 = relocated) with the following exceptions:

All live relocations are coded with 1. If a fish is relocated and is dead, then D is used. For example,

101011 - fish released on occasion 1 is relocated alive on occasions 3,5 and 6

10111D - fish released on occasion 1 is relocated alive on occasions 3,4,and 5 but is relocated dead on occasion 6.

New releases are allowed to occur on multiple occasions. The capture history of newly-released individuals should be coded with a zero (0) for the occasions before their release.

100110 - fish released on occasion 1 is relocated live on occasion 4 and 5

101000 - fish released on occasion 1 is relocated live on occasion 3

010111 - fish released on occasion 2 is relocated live on occasion 4, 5 and 6

011000 - fish released on occasion 2 is relocated live on occasion 3

001101 - fish released on occasion 3 is relocated live on occasion 4 and 6

00100D - fish released on occasion 3 is relocated dead on occasion 6.

To censor fish from the analyses, specify E after the last live encounter. For example,

10111E000 - fish released on occasion 1 is relocated alive on occasions 3,4,and 5 but is believed to have emigrated from the area by occasion 6. The capture history before the E will be used, but the fish is not included in the virtual release in occasion 6.

All life histories are summarized to reduced m-arrays (Burnham et al. (1987: page 47, Table 1.15).

The function optim is used to find F, M and P parameters that minimize the negative log-likelihood. Only cells specified in whichlivecells and whichdeadcells are used in parameter estimation.

The logit transformation is used in the estimation process to constrain values between 0 and 1. Logit-scale estimated parameters are used to calculate Sf=1/(1+exp(-B)), Sm=1/(1+exp(-C)) and P=1/(1+exp-(D)). F and M are obtained by -log(Sf) and -log(Sm).

The standard error of Sfs, Sm, P, F and M are obtained by the delta method:

SE(Sf)=sqrt((var(B)*exp(2*B))/(1+exp(B))^4),

SE(Sm)=sqrt((var(C)*exp(2*C))/(1+exp(C))^4),

SE(P)=sqrt((var(D)*exp(2*D))/(1+exp(D))^4),

SE(F)=sqrt(SE(Sf)^2/Sf^2),

SE(M)=sqrt(SE(Sm)^2/Sm^2).

All summary statistics follow Burnham and Anderson (2002). Model degrees of freedom are calculated as nlive+ndead+nnever-nreleases-1-npar where nlive is the number of whichlivecells cells, ndead is the number of whichdeadcells cells, nnever is the number of never-seen cells, nreleases is the number of releases and npar is the number of estimated parameters. Total chi-square is calculated by summing the cell chi-square values.

The program allows the configuration of different model structures (biological realistic models) for the estimation of fishing and natural mortalities, and detection probabilities. These structures are specified in Fdesign, Mdesign and Pdesign. Consider the following examples:

Example 1

Tags are relocated over seven occasions. One model structure might be constant fishing mortality estimates over occasions 1-3 and 4-6, one constant estimate of natural mortality for the entire sampling period, and one estimate of probability of detection for each occasion. To specify this model structure: Fdesign is c(“1”,“4”), Mdesign is c(“1”) and the Pdesign is c(“2:2”).

Note: The structures of Fdesign and Mdesign must always start with the first occasion, whereas the structure for Pdesign must always start with the second occasion.

Use the multiplication sign to specify occasions whose estimates of F, M or P will be taken from values of other occasions.

Example 2

Tags are relocated over six occasions. One model structure might be separate fishing mortality estimates for occasions 1-3 but assign the same parameter estimates to occasions 4-6, one constant estimate of natural mortality for occasions 1-5 and 6, and one constant probability of detection over all occasions. The Fdesign is c(“1:3*4:6”), the Mdesign is c(“1”,“6”) and the Pdesign is c(“2”).

Example 3

Specification of model 18 listed in Table 1 of Hightower et al. (2001) is shown. Each occasion represented a quarter of the year. The quarter-year design for F, M and P specifies that quarterly estimates are made in each year. Fdesign is c(“1”,“4”,“7”,“11”,“14”,“17”, “20”,“24”). Mdesign is c(“1”,“4”,“7”,“11”,“14”,“17”,“20”,“24”) and the Pdesign is c(“2”,“4”,“7”,“11”,“14”,“17”, “20”, “24”).

If the number of occasions to be assigned parameters from other occasions are less than the number of original parameters (e.g., c(“11:13*24:25”), then only the beginning sequence of original parameters equal to the number of occasions are used. For instance, in c(“11:13*24:25”), only parameters 11 and 12 would be assigned to occasions 24 and 25.

If the number of occasions to be assigned parameters from other occasions are greater than the number of original parameters (e.g., c(“11:12*24:26”)), then the last original parameter is re-cycled. In the example c(“11:12*24:26”), the parameter for occasion 12 is assigned to occasions 25 and 26.

To assist with the parameter structures, function fm_checkdesign may be used to check the desired design before use in this function.

If values of standard error are NA in the output, the hessian matrix used to claculate the variance-covariance matrix could not be inverted. If this occurs, try adjusting the reltol argument (for more information, see function optim).

In this function, the never-seen expected number is calculated by summing the live and dead probabilities, subtracting the number from 1, and then multiplying it by the number of releases. No rounding occurs in this function.

The multinomial likelihood includes the binomial coefficient.

Model averaging of model can be accomplished using the function fm_model_avg.

Note: In Hightower et al.'s original analysis, the cell probability code in SURVIV for the dead relocation in release occasion 6 had an error. The corrected analysis changed the estimates for occasions 11-13 compared to the original published values.

Value

List containing summary statistics for the model fit, model convergence status, parameter estimates estimates of fishing mortality, natural mortality, and probabilties of detection and standard errors by occasion, the parameter structure (Fdeisgn, Mdesign and Pdesign), the m-arrays, the expected (predicted) number of live and dead relocations, cell chi-square and Pearson values for live and dead relocations, matrices with the probability of being relocated alive and dead by occasion, the whichlivecells and whichdeadcells structures, and configuration label (type) used in the fm_model_avg function.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Burnham, K. P. and D. R. Anderson. 2002. Model selection and multimodel inference : A Practical Information-Theorectic Approach, 2nd edition. Spriner-Verlag, New York, NY. 488 p.

Burnham, K. P. D. R. Anderson, G. C. White, C. Brownie, and K. H. Pollock. 1987. Design and analysis methods for fish survival experiments based on release-recapture. American FIsheries Society Monograph 5, Bethesda, Maryland.

Hightower, J. E., J. R. Jackson, and K. H. Pollock. 2001. Use of telemetry methods to estimate natural and fishing mortality of striped bass in Lake Gaston, North Carolina. Transactions of the American Fisheries Society 130: 557-567.

See Also

fm_model_avg,fm_checkdesign

Examples

## Not run: 
# Set up for Full model of Hightower et al.(2001)
data(Hightower)
fm_telemetry(filetype=1,caphistory=Hightower$caphistory, Fdesign=c("1:26"),
 Mdesign=c("1:26"),Pdesign = c("2:25"),
whichlivecells=list(c(1,5,4), c(6,6,5),
 c(7,26,4)), 
whichdeadcells=list(c(1,5,4), c(6,6,6),
 c(7,26,4)),
initial=c(0.05,0.02,0.8),
control=list(reltol=1e-5,maxit=1000000))

#Set up for best model F(Qtr,yr), M constant, Pocc
fm_telemetry(filetype=1,caphistory=Hightower$caphistory, Fdesign=c("1", "4", "7", "11",
 "14", "17", "20", "24"), 
Mdesign=c("1"), Pdesign = c("2:27"),
whichlivecells=list(c(1,5,4), c(6,6,5),
 c(7,26,4)),
whichdeadcells=list(c(1,5,4), c(6,6,6),
 c(7,26,4)), 
initial=c(0.05,0.02,0.8),
 control=list(reltol=1e-8,maxit=1000000))

## End(Not run)

Description

Calculates fishing power correction ratios between two vessels or gears

Usage

fpc(cpue1 = NULL, cpue2 = NULL, method = c(1,2,3,4),  deletezerosets = FALSE, 
kapp_zeros = "paired", boot_type = "paired", nboot = 1000, dint = c(1e-9,5),
 rint = c(1e-9, 20), decimals = 2, alpha = 0.05)

Arguments

Details

The four methods for estimating fishing power correction factors given in Wilderbuer et al. (1998) are encoded.

If paired CPUE observations are both zero, the row is automatically eliminated. If deletezerosets = TRUE, the paired CPUE observations with any CPUE = 0 will be eliminated.

Zeroes are allowed in methods 1, 2 and 3.

For the Kappenman method (method=4), only non-zero CPUEs are allowed. Use kapp_zeros to select the elimination method. An unequal number of observations between vessels is allowed in this method and can result using kapp_zeros = "ind". FPC is derived by using the methodology where r that minimizes the sum of squares under the first conjecture relative to the second is estimated (Kappenman 1992: 2989; von Szalay and Brown 2001).

Standard errors and confidence intervals of FPC estimates are derived for most methods by using an approximation formula (where applicable), jackknifing and/or bootstrapping. Specify the type of bootstrapping through boot_type. For methods 1-3, jackknife estimates are provided only when boot_type="paired". If method = 4, jackknife estimates are provided only when boot_type="paired" and kapp_zeros="paired".

Confidence intervals are provided for the approximation formulae specified in Wilderbuer et al (1998), the jackknife estimates and bootstrap estimates. Confidence intervals for the jackknife method are calculated using the standard formula (estimate+/-z[alpha/2]*jackknife standard error). Bootstrap confidence intervals are derived using the percentile method (Haddon 2001).

Value

A dataframe containing method name, sample size for cpue1 (n1) and cpue2 (n2) ,mean cpue1, mean cpue2, fishing power correction (FPC), standard error from approximation formulae (U_SE), standard error from jackknifing (Jack_SE), standard error from bootstrapping (Boot_SE), lower and upper confidence intervals from approximation formulae (U_X%_LCI and U_X%_UCI),lower and upper confidence intervals from jackknifing (Jack_X%_LCI and Jack_X%_UCI) and lower and upper confidence intervals from bootstrapping (Boot_X%_LCI and Boot_X%_UCI).

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Haddon, M. 2001. Modelling and Quantitative Methods in Fisheries. Chapman & Hall/CRC Press. Boca Raton, Florida.

Kappenman, R. F. 1992. Robust estimation of the ratio of scale parameters for positive random variables. Communications in Statistics, Theory and Methods 21: 2983-2996.

von Szalay, P. G. and E. Brown. 2001. Trawl comparisons of fishing power differences and their applicability to National Marine Fisheries Service and Alask Department of Fish and Game trawl survey gear. Alaska Fishery Research Bulletin 8(2):85-95.

Wilderbuer, T. K., R. F. Kappenman and D. R. Gunderson. 1998. Analysis of fishing power correction factor estimates from a trawl comparison experiment. North American Journal of Fisheries Management 18:11-18.

Examples

## Not run: 
 #FPC for flathead sole from von Szalay and Brown 2001
   data(sole)
   fpc(cpue1=sole$nmfs,cpue2=sole$adfg,boot_type="unpaired",kapp_zeros="ind",method=c(4),
            alpha=0.05)
## End(Not run)

Description

This function finds unusual spaces or gaps in a vector of random samples

Usage

gap(x = NULL)

Arguments

Details

Values (x) are sorted from smallest to largest. Then Z values are calculated as follows: Zn-i+1=[i*(n-i)(Xn-i+1 - Xn-i)]^0.5

where n is the sample size

for i = 2,...,n calulate the 25 percent trimmed mean and divide into Z. This standardizes the distribution of the weighted gaps around a middle value of one. Suspiciously large observations should correspond to large standardized weighted gaps.

Value

vector of standardized weighted gaps

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Tukey, J. W. 1971. Exploratory data analysis. Addison-Wesley, Reading, MA. 431 pp.

Examples

y<-c(rnorm(10,10,2),1000)
 gap(y)

Description

The Gerking data frame has 14 rows and 3 columns. Marked and released sunfish in an Indiana lake for 14 days by Gerking (1953) as reported by Krebs (1989, Table 2.1).

Usage

Gerking

Format

This data frame contains the following columns:

C

column of number of captures (column names is unnecessary).

R

column of number of recaptures (column name is unnecessary).

nM

column of number of newly marked animal (column name is unnecessary).

Source

Krebs, C. J. 1989. Ecological Methodologies. Harper and Row, New York, NY. 654 p.

Description

The goosefish data frame has 40 rows and 3 columns. The mean lengths (mlen) by year and number (ss) of observations for length>=smallest length at first capture (Lc) for northern goosefish used in Gedamke and Hoenig (2006)

Usage

goosefish

Format

This data frame contains the following columns:

year

year code

mlen

mean length of goosefish, total length (cm)

ss

number of samples used to calculate mean length

Source

Gedamke, T. and J. M. Hoenig. 2006. Estimating mortality from mean length data in nonequilibrium situations, with application to the assessment of goosefish. Trans. Am. Fish. Soc. 135:476-487

Description

This function estimates parameters of Francis (1988)'s growth model using tagging data. The data are fitted using a constrained maximum likelihood optimization performed by optim using the "L-BFGS-B" method.

Usage

grotag(L1 = NULL, L2 = NULL, T1 = NULL, T2 = NULL, alpha = NULL, beta = NULL, 
 design = list(nu = 0, m = 0, p = 0, sea = 0), 
 stvalue = list(sigma = 0.9, nu = 0.4, m = -1, p = 0.1, u = 0.4, w = 0.4), 
 upper = list(sigma = 5, nu = 1, m = 2, p = 1, u = 1, w = 1), 
 lower = list(sigma = 0, nu = 0, m = -2, p = 0, u = 0, w = 0), gestimate = TRUE, 
 st.ga = NULL, st.gb = NULL, st.galow = NULL, st.gaup = NULL, st.gblow = NULL,
 st.gbup = NULL, control = list(maxit = 10000))

Arguments

Details

The methods of Francis (1988) are used on tagging data to the estimate of growth and growth variability. The estimation of all models discussed is allowed. The growth variability defined by equation 5 in the reference is used throughout.

Value

Author(s)

Marco Kienzle [email protected]

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

Francis, R.I.C.C., 1988. Maximum likelihood estimation of growth and growth variability from tagging data. New Zealand Journal of Marine and Freshwater Research, 22, p.42-51.

See Also

grotagplus

Examples

data(bonito)

#Model 4 of Francis (1988)
with(bonito,
 grotag(L1=L1, L2=L2, T1=T1, T2=T2,alpha=35,beta=55,
 	design=list(nu=1,m=1,p=1,sea=1),
 	stvalue=list(sigma=0.9,nu=0.4,m=-1,p=0.2,u=0.4,w=0.4),
 	upper=list(sigma=5,nu=1,m=2,p=0.5,u=1,w=1),
 	lower=list(sigma=0,nu=0,m=-2,p=0.0,u=0,w=0),control=list(maxit=1e4)))

Description

This is an extension of fishmethods function grotag to allow a wider variety of growth models and also the simultaneous analysis of multiple tagging datasets with parameter sharing between datasets (see Details).

As in grotag, the data are fitted using a constrained maximum likelihood optimization performed by optim using the "L-BFGS-B" method. Estimated parameters can include galpha, gbeta (mean annual growth at reference lengths alpha and beta); b (a curvature parameter for the Schnute models); Lstar (a transitional length for the asymptotic model); m, s (mean and s.d. of the measurement error for length increment); nu, t (growth variability); p (outlier probability); u, w (magnitude and phase of seasonal growth).

Usage

grotagplus(tagdata, dataID=NULL,alpha, beta = NULL,
 model=list(mean="Francis",var="linear",seas="sinusoid"),
 design, stvalue, upper, lower,fixvalue=NULL,
 traj.Linit=c(alpha,beta),control = list(maxit = 10000), debug = FALSE)

Arguments

Details

Valid values of model$mean are "Francis" as in Francis (1988). "Schnute" as in Francis (1995). "Schnute.aeq0" special case of Schnute - see equns (5.3), (5.4) of Francis (1995). "asymptotic" as in Cranfield et al. (1996).

Valid values of model$var are "linear" as used in the example in Francis(1988) - see equn (5). "capped" as in equn (6) of Francis(1988). "exponential" as in equn (7) of Francis(1988).

"asymptotic" as in equn (8) of Francis(1988). "least-squares" ignore individual variability and fit data by least-squares, as in Model 1 of Francis(1988).

Valid values of model$seas are "sinusoid" as in model 4 of Francis(1988). "switched" as in Francis & Winstanley (1989). "none" as in all but model 4 of Francis(1988).

The option of multiple data sets with parameter sharing is intended to allow for the situation where we wish to estimate different mean growth for two or more datasets but can reasonably assume that other parameters (e.g., for growth variability, measurement error, outlier contamination) are the same for all datasets. This should produces stronger estimates of these other parameters. For example, Francis & Francis (1992) allow growth to differ by sex, and in Francis & Winstanley (1989) it differs by stock and/or habitat.

grotagplus may fail if parameter starting values are too distant from their true value, or if parameter bounds are too wide. Try changing these values. Sometimes reasonable starting values can be found by fitting the model with other parameters fixed at plausible values.

Value

Author(s)

Chris Francis [email protected]

Marco Kienzle [email protected]

Gary A. Nelson, Massachusetts Division of Marine Fisheries [email protected]

References

1 Francis, R.I.C.C., 1988. Maximum likelihood estimation of growth and growth variability from tagging data. New Zealand Journal of Marine and Freshwater Research, 22, p.42-51.

2 Cranfield, H.J., Michael, K.P., and Francis, R.I.C.C. 1996. Growth rates of five species of subtidal clam on a beach in the South Island, New Zealand. Marine and Freshwater Research 47: 773-784.

3 Francis, R.I.C.C. 1995. An alternative mark-recapture analogue of Schnute"s growth model. Fisheries Research 23: 95-111.

4 Francis, R.I.C.C. and Winstanley, R.H. 1989. Differences in growth rates between habitats of southeast Australian snapper (Chrysophrys auratus). Australian Journal of Marine & Freshwater Research 40: 703-710.

5 Francis, M.P. and Francis, R.I.C.C. 1992. Growth rate estimates for New Zealand rig (Mustelus lenticulatus). Australian Journal of Marine and Freshwater Research 43: 1157-1176.

See Also

plot.grotagplus print.grotagplus

Examples

#Model 4 of Francis (1988)
data(bonito)
grotagplus(bonito,alpha=35,beta=55,
               design=list(galpha=1,gbeta=1,s=1,nu=1,m=1,p=1,u=1,w=1),
               stvalue=list(s=0.81,nu=0.3,m=0,p=0.01,u=0.5,w=0.5),
               upper=list(s=3,nu=1,m=2,p=0.1,u=1,w=1),
               lower=list(s=0.1,nu=0.1,m=-2,p=0,u=0,w=0))

#Model 1 of Francis (1988), using least-squares fit
grotagplus(bonito,alpha=35,beta=55,
               model=list(mean="Francis",var="least-squares",seas="none"),
               design=list(galpha=1,gbeta=1,s=1,p=0),
               stvalue=list(s=1.8),upper=list(s=3),lower=list(s=1))

#Paphies donacina model in Table 4 of Cranfield et al (1996) with
#asymptotic model
data(P.donacina)
grotagplus(P.donacina,alpha=50,beta=80,
       model=list(mean="asymptotic",var="linear",seas="none"),
       design=list(galpha=1,gbeta=1,Lstar=0,s=1,nu=0,m=0,p=0),
       stvalue=list(galpha=10,gbeta=1.5,s=2),
       upper=list(galpha=15,gbeta=2.7,s=4),
       lower=list(galpha=7,gbeta=0.2,s=0.5),
       fixvalue=list(Lstar=80))

#Paphies donacina model in Table 4 of Cranfield et al (1996) with
#asymptotic model
data(P.donacina)
grotagplus(P.donacina,alpha=50,beta=80,
       model=list(mean="asymptotic",var="linear",seas="none"),
       design=list(galpha=1,gbeta=1,Lstar=0,s=1,nu=0,m=0,p=0),
       stvalue=list(galpha=10,gbeta=1.5,s=2),
       upper=list(galpha=15,gbeta=2.7,s=4),
       lower=list(galpha=7,gbeta=0.2,s=0.5),
       fixvalue=list(Lstar=80))

# Model 4 fit from Francis and Francis (1992) with different growth by sex
data(rig)
grotagplus(rig,dataID="Sex",alpha=70,beta=100,
           model=list(mean="Francis",var="linear",seas="none"),
          design=list(galpha=list("F","M"),gbeta=list("F","M"),s=1,nu=1,m=0,p=0),
          stvalue=list(galpha=c(5,4),gbeta=c(3,2),s=2,nu=0.5),
          upper=list(galpha=c(8,6),gbeta=c(5,4),s=4,nu=1),
          lower=list(galpha=c(3,2),gbeta=c(1.5,1),s=0.5,nu=0.2))

#Example where all parameters are fixed 
# to the values estimated values for model 4 of Francis and Francis (1992)]
grotagplus(rig,dataID="Sex",alpha=70,beta=100,
          model=list(mean="Francis",var="linear",seas="none"),
          design=list(galpha=0,gbeta=0,s=0,nu=0,m=0,p=0),
          stvalue=list(),upper=list(),lower=list(),
          fixvalue=list(galpha=list(design=list("F","M"),value=c(5.87,3.67)),
          gbeta=list(design=list("F","M"),value=c(2.52,1.73)),s=1.57,nu=0.58))