Publication - Progress report

Scottish Marine and Freshwater Science Vol 6 No 15: Spatial dynamics of scallops in relation to the Orkney dive fishery. Report of Fishing Industry Science Alliance (FISA) Project 03/12

Published: 25 November 2015
Topic: Marine and fisheries
ISBN: 9781785448294

This report details the results of studies which have been done under Fishing Industry Science Alliance (FISA) project 03/12 which quantify the rate of spatial turnover in a scallop population at a small spatial scale in Orkney, and understanding the exte

Supporting documents

Appendix 1

A Simple Maximum Likelihood Model for Depletion Fishing

For repeated passes over the same plot the expected catch number on each pass Ĉ _i is a function of initial population size N and capture probabilities P _i on each occasion ( i=1 to np):

Equation

Assuming a Poisson distribution, the total log-likelihood for catches over the np passes is:

Equation

This model has been implemented in the R statistical package (R Core Team, 2014) using the 'maxLik' library (Henningson & Toomet, 2011) for maximizing this log-likelihood and the 'numDeriv' library (Gilbert & Varadhan, 2012) to estimate numerical gradients and Hessian matrix. Initial population size and capture probabilities are constrained within feasible bounds by logarithmic and logistic transforms respectively. The next page shows an R script for fitting the model to the Wyre depletion data, with capture probability assumed to remain constant between passes. This code allows simultaneous estimation for multiple groups, with between-group constraints on capture probabilities defined through a GLM-like design matrix. Five models are shown defined in the code below:

(i) unconstrained estimation of model parameters separately for each group;

(ii) capture probability constrained to be an additive logistic function of size-group and fishing occasion effects;

(iii) capture probability constrained to be the same in each size-group but differing between occasions;

(iv) capture probability constrained to be the same on each occasion, but differing between size-groups; and

(v) capture probability constrained to be constant across occasions and size-groups.

Gilbert, P. & Varadhan, R., 2012. numDeriv: Accurate Numerical Derivatives. R package version 2012.9-1.

Henningsen, A. & Toomet, O., 2011. maxLik: A package for maximum likelihood estimation in R. Computational Statistics, 26, 443-458.

R Core Team, 2014. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

R Script for Defining and Fitting Maximum Likelihood Depletion Models

##----------------------------------------------------------------

## MAXIMUM-LIKELIHOOD DEPLETION MODEL

##----------------------------------------------------------------

require(MASS) ## for generalized inverse function, ginv()

require(Matrix) ## for rankMatrix() function

require(maxLik) ## maximum likelihood estimation

require(numDeriv) ## numerical derivatives for gradient & Hessian

## N.B. numDeriv must be loaded AFTER maxLik so that hessian

## function from maxLik is masked by that from numDeriv rather

## than vice versa

##----------------------------------------------------------------

## POISSON LOG-LIKELIHOOD FUNCTION FOR A DEPLETION MODEL

## p[npp+ng] is the estimated parameter vector

## C[nr] is the catch vector

## Pdesign[nr,npp] is the design matrix for P parameters

## gt[nr,2] gives group and time qualifiers

## npp is the number of P parameters in the estimated vector

## nr is the number of data rows

## ng is the number of groups (populations)

## ntmax is the maximum number of passes for any population

##----------------------------------------------------------------

loglikDeplete<-function(p,C,Pdesign,gt,npp,nr,ng,ntmax)

{

##Extract P from estimated parameters

P<-matrix(nrow=ng,ncol=ntmax)

for(i in 1:nr) {

g<-gt[i,1]

t<-gt[i,2]

P[g,t]<-0

for(j in 1:npp) {

P[g,t]<-P[g,t]+Pdesign[i,j]*p[j]

}

P[g,t]<-1/(1+exp(-P[g,t]))

}

## Extract N from estimated parameters

N<-vector(length=ng)

for(i in 1:ng) {

N[i]<-exp(p[npp+i])

}

## form log-likelihood

LL<-0

for(i in 1:nr) {

g<-gt[i,1]

t<-gt[i,2]

nj<-t-1

Pprod<-1

if(nj>0){

for(j in 1:nj) {

Pprod<-Pprod*(1-P[g,j])

}

Cfit<-N[g]*Pprod*P[g,t]

LL<-LL+C[i]*log(Cfit)-Cfit

}

return(LL)

}

##----------------------------------------------------------------

## CALCULATE FITTED VALUES

## Same calculations as loglikDeplete, but returns fitted

## values of catch numbers rather than log-likelihood

##----------------------------------------------------------------

calcFit<-function(p,C,Pdesign,gt,npp,nr,ng,ntmax)

{

##Extract P from estimated parameters

P<-matrix(nrow=ng,ncol=ntmax)

for(i in 1:nr) {

g<-gt[i,1]

t<-gt[i,2]

P[g,t]<-0

for(j in 1:npp) {

P[g,t]<-P[g,t]+Pdesign[i,j]*p[j]

}

P[g,t]<-1/(1+exp(-P[g,t]))

}

## Extract N from estimated parameters

N<-vector(length=ng)

for(i in 1:ng) {

N[i]<-exp(p[npp+i])

}

## form log-likelihood

Cfit<-vector(length=nr)

for(i in 1:nr) {

g<-gt[i,1]

t<-gt[i,2]

nj<-t-1

Pprod<-1

if(nj>0){

for(j in 1:nj) {

Pprod<-Pprod*(1-P[g,j])

}

Cfit[i]<-N[g]*Pprod*P[g,t]

}

return(Cfit)

}

##----------------------------------------------------------------

## NUMERICAL GRADIENT OF THE LOG-LIKELIHOOD

##----------------------------------------------------------------

gradDeplete<-function(p,C,Pdesign,gt,npp,nr,ng,ntmax)

{

g<-grad(loglikDeplete,p,C=C,Pdesign=Pdesign,gt=gt,npp=npp,

nr=nr,ng=ng,ntmax=ntmax)

return(g)

}

##----------------------------------------------------------------

## NUMERICAL HESSIAN OF THE LOG-LIKELIHOOD

##----------------------------------------------------------------

hessDeplete<-function(p,C,Pdesign,gt,npp,nr,ng,ntmax)

{

h<-hessian(loglikDeplete,p,C=C,Pdesign=Pdesign,gt=gt,npp=npp,

nr=nr,ng=ng,ntmax=ntmax)

return(h)

}

##----------------------------------------------------------------

## FIT THE DEPLETION MODEL BY MAXIMUM LIKELIHOOD

## Returns dmod as object of class maxLike (estimated parameters

## are returned as component dmod$estimate)

## Returns var as covariance matrix

## Returns se as vector of standard errors for model parameters

## Returns NP as number of identifiable model parameters

## Returns AIC as Akaike Information Criterion

## Returns AICc as small sample AIC

##-----------------------------------------------------------------

FitDeplete<-function(p,C,Pdesign,gt,npp,nr,ng,ntmax)

{

## Fit the model

dmod<-maxLik(loglikDeplete,grad=gradDeplete,hess=hessDeplete,

start=p,method="BFGS",

C=C,Pdesign=Pdesign,gt=gt,npp=npp,nr=nr,ng=ng,

ntmax=ntmax)

## Get model fit statistics

var<-ginv(-dmod$hessian)

se<-sqrt(diag(var))

NP<-rankMatrix(dmod$hessian,method="maybeGrad")

attributes(NP)<-NULL ## strip the object down to a value

AIC<-(-2*dmod$maximum)+2*NP

AICc<-AIC+2*NP*(NP+1)/(length(C)-NP-1)

## Extract logit(P) estimates and SEs

logitP<-vector(length=nr)

SElogitP<-vector(length=nr)

for(i in 1:nr) {

g<-gt[i,1]

t<-gt[i,2]

logitP[i]<-0

SElogitP[i]<-0

for(j in 1:npp) {

logitP[i]<-logitP[i]+Pdesign[i,j]*dmod$estimate[j]

for(k in 1:npp) {

SElogitP[i]<-SElogitP[i]+Pdesign[i,j]*Pdesign[i,k]*var[j,k]

}

SElogitP[i]<-sqrt(SElogitP[i])

}

## Extract log(N) estimates and SEs

logN<-vector(length=ng)

SElogN<-vector(length=ng)

for(i in 1:ng) {

logN[i]<-dmod$estimate[npp+i]

SElogN[i]<-se[npp+i]

}

## Get the model deviance

LLsaturated<-0

for(i in 1:nr) {

if(C[i]>0) {

LLsaturated<-LLsaturated+C[i]*log(C[i])-C[i]

}

deviance<-2*LLsaturated-2*dmod$maximum

## Get the model fit

Cfit<-vector(length=nr)

Cfit<-calcFit(p=dmod$estimate,C=C,Pdesign=Pdesign,gt=gt,npp=npp,

nr=nr,ng=ng,ntmax=ntmax)

return(list(dmod=dmod,var=var,se=se,logitP=logitP,

SElogitP=SElogitP,logN=logN,SElogN=SElogN,

NP=NP,AIC=AIC,AICc=AICc,deviance=deviance,

Cfit=Cfit))

}

##-----------------------------------------------------------------

## DATA FOR WYRE EXPERIMENT

ng<-15 ## 5 occasions x 3 size groups

nr<-45 ## 15 groups x 3 passes per group

ntmax<-3 ## maximum of 3 passes

## Catch vector

C<-as.vector(c(41,15,13,

56,6,0,

76,6,4,

47,32,34,

53,15,10,

8,12,0,

57,21,3,

94,27,4,

11,5,1,

32,30,25,

112,43,12,

39,1,1,

10,31,9,

43,22,12,

12,2,2),

mode="numeric")

## Group and time qualifiers

gt<-matrix(c(1,1,

1,2,

1,3,

2,1,

2,2,

2,3,

3,1,

3,2,

3,3,

4,1,

4,2,

4,3,

5,1,

5,2,

5,3,

6,1,

6,2,

6,3,

7,1,

7,2,

7,3,

8,1,

8,2,

8,3,

9,1,

9,2,

9,3,

10,1,

10,2,

10,3,

11,1,

11,2,

11,3,

12,1,

12,2,

12,3,

13,1,

13,2,

13,3,

14,1,

14,2,

14,3,

15,1,

15,2,

15,3),

nrow=nr,ncol=2,byrow=TRUE)

##-----------------------------------------------------------------

## SET UP MODEL STRUCTURES - Size Group * Occasion

npp=15 ## P parameters for 15 populations

## Initial values for parameter estimates

## first 15 are logit(P), next 15 are log(N)

p<-as.vector(c(rep(0,15),rep(4,15)))

## Design matrix mapping estimated to structural parameters

Pdesign<-matrix(c(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,

0,0,0,0,0,0,0,0,0,0,0,0,0,0,1),

nrow=nr,ncol=npp,byrow=TRUE)

##-----------------------------------------------------------------

## RUN THE MODEL - Size Group * Occasion

DM1<-FitDeplete(p,C,Pdesign,gt,npp,nr,ng,ntmax)

##-----------------------------------------------------------------

## SET UP MODEL STRUCTURES - Size Group + Occasion

npp=7 ## P parameters

## Initial values for parameter estimates

## first 7 are logit(P), next 15 are log(N)

p<-as.vector(c(rep(0,7),rep(4,15)))

## Design matrix mapping estimated to structural parameters

Pdesign<-matrix(c(1,0,0,0,0,0,0,

1,0,0,0,0,0,0,

1,0,0,0,0,1,0,

1,0,0,0,0,0,1,

0,1,0,0,0,0,0,

0,1,0,0,0,1,0,

0,1,0,0,0,0,1,

0,0,1,0,0,0,0,

0,0,1,0,0,1,0,

0,0,1,0,0,0,1,

0,0,0,1,0,0,0,

0,0,0,1,0,1,0,

0,0,0,1,0,0,1,

0,0,0,0,1,0,0,

0,0,0,0,1,1,0,

0,0,0,0,1,0,1,

0,0,0,0,1,0,1),

nrow=nr,ncol=npp,byrow=TRUE)

##-----------------------------------------------------------------

## RUN THE MODEL - Size Group + Occasion

DM2<-FitDeplete(p,C,Pdesign,gt,npp,nr,ng,ntmax)

##-----------------------------------------------------------------

## SET UP MODEL STRUCTURES - Occasion

npp=5 ## P parameters

## Initial values for parameter estimates

## first 5 are logit(P), next 15 are log(N)

p<-as.vector(c(rep(0,5),rep(4,15)))

## Design matrix mapping estimated to structural parameters

Pdesign<-matrix(c(1,0,0,0,0,

1,0,0,0,0,

0,1,0,0,0,

0,0,1,0,0,

0,0,0,1,0,

0,0,0,0,1,

0,0,0,0,1),

nrow=nr,ncol=npp,byrow=TRUE)

##-----------------------------------------------------------------

## RUN THE MODEL - Occasion

DM3<-FitDeplete(p,C,Pdesign,gt,npp,nr,ng,ntmax)

##-----------------------------------------------------------------

## SET UP MODEL STRUCTURES - Size Group

npp=3 ## P parameters

## Initial values for parameter estimates

## first 3 are logit(P), next 15 are log(N)

p<-as.vector(c(rep(0,3),rep(4,15)))

## Design matrix mapping estimated to structural parameters

Pdesign<-matrix(c(1,0,0,

1,0,0,

0,1,0,

0,0,1,

1,0,0,

0,1,0,

0,0,1,

1,0,0,

0,1,0,

0,0,1,

1,0,0,

0,1,0,

0,0,1,

1,0,0,

0,1,0,

0,0,1,

0,0,1),

nrow=nr,ncol=npp,byrow=TRUE)

##-----------------------------------------------------------------

## RUN THE MODEL - Size Group

DM4<-FitDeplete(p,C,Pdesign,gt,npp,nr,ng,ntmax)

##-----------------------------------------------------------------

## SET UP MODEL STRUCTURES - Constant P

npp=1 ## P parameters

## Initial values for parameter estimates

## first 1 is logit(P), next 15 are log(N)

p<-as.vector(c(0,rep(4,15)))

## Design matrix mapping estimated to structural parameters

Pdesign<-matrix(c(rep(1,45)),

nrow=nr,ncol=npp,byrow=TRUE)

##-----------------------------------------------------------------

## RUN THE MODEL - Constant P

DM5<-FitDeplete(p,C,Pdesign,gt,npp,nr,ng,ntmax)

Contact

Stephanie Dickens

There is a problem

Thanks for your feedback

Was this helpful?

Yes

Your comments

Your feedback helps us to improve this website. Do not give any personal information because we cannot reply to you directly.

Choose a reason for your feedback

Your comments

Your feedback helps us to improve this website. Do not give any personal information because we cannot reply to you directly.

Yes, but

Choose a reason for your feedback

Your comments

Your feedback helps us to improve this website. Do not give any personal information because we cannot reply to you directly.