Package 'MQMF'

Description

'%ni%' identifies which element in x is NOT in y

Usage

x %ni% y

Arguments

Examples

x <- 1:10
   y <- 6:15
   x
   y
   x[(x %ni% y)]   # are not in y

Description

A dataset of fishery data for blacklip abalone (Haliotis rubra) from part of the Tasmanian west coast for the years 1985 - 2008. It contains a data.frame containing the year, the catch, and the standardized CPUE from four statistical blocks of Tasmania's west coast combined. In particular, it can be used when fitting a surplus production model. Workable initial parameter values, before log-transformation might be: r= 0.4, K=9400, Binit=3400, sigma=0.05 for the Schaefer version, while these also work for the Fox model one could more efficiently use r=0.3, K=12000, Binit=4000, sigma=0.05.

Format

A data.frame of three columns

year

the annual year in which the catches occurred

catch

the reported landed catch in tonnes, to the nearest kilogram

cpue

the standardized catch-per-unit-effort for this dive fishery

Subjects

• Surplus Production Modelling, Schaefer and Fox models

• Model fitting using maximum likelihood

• Uncertainty examples

Source

Catch data from Mundy, C. and J. McAllister (2019) Tasmanian abalone fishery assessment 2018, Institute for Marine and Antarctic Studies, University of Tasmania, 190p. ISBN: 978-1-925646-46-7. The cpue data is an unpublished early attempt at standardizing the cpue data with respect to month, block, and diver. Many more details are now included in such analyses.

Examples

data(abdat)
 print(abdat)
 oldpar <- par(no.readonly=TRUE)
 plot(abdat$year,abdat$cpue, type="l",xlab="year",ylab="CPUE",
      panel.first=grid())
 points(abdat$year,abdat$cpue,pch=16,cex=1.2)
 par(oldpar)

Description

addcontours is used to add contours to a dense plot of xy points such as might be generated when conducting an analysis of the the uncertainty associated with a stock assessment, or other analysis using a bootstrap, a Bayesian MCMC, or even using asymptotic errors and sampling from a multi-variate normal distribution. addcontours first uses the kde2d function from the MASS package to translate the density of points into 2-D kernal densities, and then searches through the resulting densities for those points that would identify approximate contours. Finally it calls the contour function to add the identified approximate contours to the xy plot.

Usage

addcontours(
  xval,
  yval,
  xrange,
  yrange,
  ngrid = 100,
  contval = c(0.5, 0.9),
  lwd = 1,
  col = 1
)

Arguments

Value

nothing but it does add contours to a plot of points

Examples

library(mvtnorm)
 library(MASS)
 data(abdat)
 param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05)) 
 bestmod <- nlm(f=negLL,p=param,funk=simpspm,indat=abdat, 
                hessian=TRUE,logobs=log(abdat$cpue),
                typsize=magnitude(param),iterlim=1000)
 optpar <- bestmod$estimate
 vcov <- solve(bestmod$hessian)      # solve inverts matrices
 columns <- c("r","K","Binit","sigma")
 N <- 1000  # the contours improve as N increases; try 5000
 mvnpar <- matrix(exp(rmvnorm(N,mean=optpar,sigma=vcov)),
                  nrow=N,ncol=4,dimnames=list(1:N,columns))
 xv <- mvnpar[,"K"]
 yv <- mvnpar[,"r"]
 # plotprep(width=6,height=5,newdev = FALSE)
 plot(xv,yv,type="p") # use default 0.5 and 0.9 contours.
 addcontours(xv,yv,range(xv),range(yv),lwd=2,col=2)
 points(mean(xv),mean(yv),pch=16,cex=1.5,col=2)

Description

addlnorm estimates a log-normal distribution from the output of a histogram of a data set.

Usage

addlnorm(inhist, xdata, inc = 0.01)

Arguments

Value

a 4 x N matrix of x and y values to be used to plot the fitted normal probability density function.Combined with estimates of mean(log(indata)) and log(sd(indata))

Examples

oldpar <- par(no.readonly=TRUE)
 egdata <- rlnorm(200,meanlog=0.075,sdlog=0.5)
 outh <- hist(egdata,main="",col=2,breaks=seq(0,8,0.2))
 ans <- addlnorm(outh,egdata)
 lines(ans[,"x"],ans[,"y"],lwd=2,col=4) 
 par(oldpar)

Description

addnorm adds a normal distribution to a histogram of a data set. This is generally to be used to illustrate whether log-transformation normalizes a set of catch or cpue data.

Usage

addnorm(inhist, xdata, inc = 0.01)

Arguments

Value

a list with a vector of 'x' values and a vector of 'y' values (to be used to plot the fitted normal probability density function), and a vector called 'stats' containing the mean and standard deviation of the input data

Examples

oldpar <- par(no.readonly=TRUE)
 x <- rnorm(1000,mean=5,sd=1)
 #plotprep(height=6,width=4,newdev=FALSE)
 par(mfrow= c(1,1),mai=c(0.5,0.5,0.3,0.05))
 par(cex=0.75, mgp=c(1.5,0.35,0), font.axis=7)
 outH <- hist(x,breaks=25,col=3,main="")
 nline <- addnorm(outH,x)
 lines(nline$x,nline$y,lwd=3,col=2)
 print(nline$stats)
 par(oldpar)

Description

aicbic calculates and returns the AIC and BIC using the standard definitions. It defaults to assuming that negative log- likelihoods have been used in the model fitting, but provides the option of having used SSQ (set nLL to FALSE). If using SSQ it uses Burnham and Anderson's (2002) definition but sets BIC to NA. aicbic can recognize the outputs from optim, nlm, and nlminb.

Usage

aicbic(model, dat, nLL = TRUE)

Arguments

Value

a vector of four numbers, AIC first, then BIC, then negLL or SSQ, depending on nLL, then number of parameters p

References

Burnham, K.P. and D.R. Anderson (2002) Model Selection and Inference. A Practical Information-Theoretic Approach. Second Edition Springer-Verlag, New York. 488 p.

Examples

data(blackisland); bi <- blackisland
param <- c(Linf=170.0,K=0.3,sigma=4.0)
modelvb <- nlm(f=negNLL,p=param,funk=fabens,observed=bi$dl,indat=bi,
               initL="l1",delT="dt") # could have used the defaults
# Do not worry about the warning messages
aicbic(modelvb,blackisland)  # 588.3382 596.3846 291.1691   3

Description

altnegLL calculates negLLM using the simplification from Haddon (2011) using the ssq calculated within the function spm

Usage

altnegLL(inp, indat)

Arguments

Value

a single value, the negative log-likelihood

Examples

data(dataspm)
 pars <- log(c(r=0.2,K=6000,Binit=2800,sigma=0.2))
 ans <- fitSPM(pars,fish=dataspm,schaefer=TRUE,maxiter=1000)
 outfit(ans)
 altnegLL(ans$estimate,dataspm) # should be -12.12879

Description

bce the Baranov Catch Equation. The total mortality of fish in an exploited population is made up of fish being killed by fishing and others dying naturally. We use the bce to estimate the catch (those killed by fishing). The bce has value because some fish that would be expected to die naturally can be expected to be caught and killed by fishing, so estimating the catch is rather more complex than numbers of fish available times the harvest rate, Nt x Ht. If is possible to use the Baranov Catch Equation when calculating the expected catches then it is usually better to do so.

Usage

bce(M, Fat, Nt, ages)

Arguments

Value

a matrix of surviving numbers-at age, total mortality-at-age, and catch-at-age

Examples

age <- 0:25
Ft <- -log(1 - 0.2) # harvest rate of 0.2
Faa <- rep(Ft,length(age))
M <- 0.12
Nt <- 1000
bce(M,Fat=Faa,Nt,ages=age)   # should give 188.8862

Description

bh implements the Beverton-Holt stock recruitment equation R = a x B/(b + B), where R is the recruitment, a and b are the parameters and B is the spawning biomass. In this parameterization, a is the maximum recruitment level and b is the biomass required to generate 0.5 x maximum recruitment. An common alternative parameterization is to use R = S/(alpha + beta x S). Do not confuse a with alpha and b with beta!

Usage

bh(p, B)

Arguments

Value

a vector, the same length as B, of the predicted recruitment(s)

Examples

B <- 1:3000
  rec <- bh(c(1000,200),B)
  plot1(B,rec,xlab="SpB",ylab="Recruitment",lwd=2)

Description

A 108 x 4 data.frame containing dt, the time in years between tagging and recapture, l1 the shell length at tagging, and l2, the length at recapture, with the growth increment, dl as the last column. This data can be used to estimate the growth characteristics of abalone from the Black Island site, which is on the south west coast of Tasmania, Australia. The mean time interval between tagging and recapture is 1 year and 1 week, 1.02 years, which reflects the practical problems of taking a vessel around the bottom of Tasmania, where it is essential to wait on suitable weather for such sub-tidal field work.

Format

A data.frame of four columns

dt

the time between tagging and recapture, in years

l1

the shell length when tagged in mm

l2

the shell length at recapture in mm

dl

the growth increment between tagging and recapture in mm; there are zero values.

Subjects

• Tagging data

• Estimation of individual growth

• Binomial likelihoods

• Faben's version of the von Bertalanffy curve

Source

Thanks to Dr Craig Mundy and the abalone team at the Institute of Marine and Antarctic Studies, of the University of Tasmania for the use of this data.

Examples

data(blackisland)
 print(head(blackisland,20))
 oldpar <- par(no.readonly=TRUE)
 plot(blackisland$l1,blackisland$dl,type="p",pch=16,
 xlab="Initial Length mm",ylab="Growth Increment mm",
 panel.first=grid())
 abline(h=0)
 par(oldpar)

Description

bracket scans through a series of predicted values for the location of a target value of the y-axis and returns the two y values that bracket the target, perhaps finding the values closest to 0.5 in a vector between 0 and 1. It also returns the x-axis values that gave rise to the two values bracketing the target, and finally returns the target. For example, imagine generating the proportion mature for a given length of fish using an equation for which there was no analytical solution to what the value of the L50 or the inter-quartile distance was. Bracket can find the two lengths that generate proportions just below 0.5 and just above. It does not matter if, by chance, the target is one of those y-axis values.

Usage

bracket(x, yaxis, xaxis)

Arguments

Value

a vector of 5 values, left, right, bottom, top and target

See Also

linter

Examples

L = seq(60,160,1)
 p=c(a=0.075,b=0.075,c=1.0,alpha=100)
 asym <- srug(p=p,sizeage=L) # Schnute and Richards unified growth curve
 L25 <- linter(bracket(0.25,asym,L)) 
 L50 <- linter(bracket(0.5,asym,L)) 
 L75 <- linter(bracket(0.75,asym,L)) 
 ans <- c(L25,L50,L75,L50-L25,L75-L50)
 {cat("   L25    L50      L75   L50-L25 L75-L50 \n")
 cat(round(ans,4),"\n")}

Description

calcprior is used to include a prior probability into Bayesian calculations. calcprior is a template for generating such priors. The default given here is to return a constant small number for the prior probability, it needs to sum to 1.0 across the replicates returned by do_MCMC. If non-uniform priors are required write a different function and in do_MCMC point priorcalc at it. Whatever function you define needs to have the same input parameters as this calcprior, i.e. the parameters and N. If something else if required then do_MCMC will need modification in the two places where priorcalc is used. Alternatively, the ellipsis, ..., might be used.

Usage

calcprior(pars, N)

Arguments

Value

the sum of a vector of small constant values to act as priors.

Examples

param <- log(c(0.4,9400,3400,0.05))  
calcprior(pars=param,N=20000)  # should give -39.61395

Description

chapter2 contains no active function but rather acts as a repository for the various example code chunks found in chapter2. There are 15 r-code chunks in chapter2.

Usage

chapter2(verbose = TRUE)

Arguments

Examples

## Not run: 
# All the example code from  # A Non-Introduction to R
### Using Functions

# R-chunk 1  Page 18
#make a function called countones2, don't overwrite original

countones2 <- function(x) return(length(which(x == 1)))  # or
countones3 <- function(x) return(length(x[x == 1]))
vect <- c(1,2,3,1,2,3,1,2,3)  # there are three ones
countones2(vect)  # should both give the answer: 3
countones3(vect)
set.seed(7100809) # if repeatability is desirable.
matdat <- matrix(trunc(runif(40)*10),nrow=5,ncol=8)
matdat #a five by eight matrix of random numbers between 0 - 9
apply(matdat,2,countones3)  # apply countones3 to 8 columns
apply(matdat,1,countones3)  # apply countones3 to 5 rows


# R-chunk 2 Page 19
 #A more complex function prepares to plot a single base graphic
 #It has syntax for opening a window outside of Rstudio and
 #defining a base graphic. It includes oldpar <- par(no.readonly=TRUE)
 #which is returned invisibly so that the original 'par' settings
 #can be recovered using par(oldpar) after completion of your plot.

plotprep2 <- function(plots=c(1,1),width=6, height=3.75,usefont=7,
                      newdev=TRUE) {
  if ((names(dev.cur()) %in% c("null device","RStudioGD")) &
      (newdev)) {
    dev.new(width=width,height=height,noRStudioGD = TRUE)
  }
  oldpar <- par(no.readonly=TRUE)  # not in the book's example
  par(mfrow=plots,mai=c(0.45,0.45,0.1,0.05),oma=c(0,0,0,0))
  par(cex=0.75,mgp=c(1.35,0.35,0),font.axis=usefont,font=usefont,
      font.lab=usefont)
  return(invisible(oldpar))
}  #  see ?plotprep; see also parsyn() and parset()


### Random Number Generation
# R-chunk 3 pages 20 - 21
#Examine the use of random seeds.

seed <- getseed()  # you will very likely get different answers
set.seed(seed)
round(rnorm(5),5)
set.seed(123456)
round(rnorm(5),5)
set.seed(seed)
round(rnorm(5),5)

### Plotting in R
# R-chunk 4  page 22
#library(MQMF)   # The development of a simple graph  see Fig. 2.1
# The statements below open the RStudio graphics window, but opening
# a separate graphics window using plotprep is sometimes clearer.

data("LatA")  #LatA = length at age data; try properties(LatA)
# plotprep(width=6.0,height=5.0,newdev=FALSE) #unhash for external plot
oldpar <- par(no.readonly=TRUE)  # not in the book's example
setpalette("R4") #a more balanced, default palette see its help
par(mfrow=c(2,2),mai=c(0.45,0.45,0.1,0.05))  # see ?parsyn
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
hist(LatA$age) #examine effect of different input parameters
hist(LatA$age,breaks=20,col=3,main="") # 3=green #try ?hist
hist(LatA$age,breaks=30,main="",col=4) # 4=blue
hist(LatA$age, breaks=30,col=2, main="", xlim=c(0,43), #2=red
     xlab="Age (years)",ylab="Count")
par(oldpar)  # not in the book's example

### Dealing with Factors
# R-chunk 5  pages 23 - 24
#Dealing with factors/categories can be tricky

DepCat <- as.factor(rep(seq(300,600,50),2)); DepCat
try(5 * DepCat[3], silent=FALSE) #only returns NA and a warning!
as.numeric(DepCat) # returns the levels not the original values
as.numeric(levels(DepCat)) #converts 7 levels not the replicates
DepCat <- as.numeric(levels(DepCat))[DepCat] # try ?facttonum
#converts replicates in DepCat to numbers, not just the levels
5 * DepCat[3]   # now treat DepCat as numeric
DepCat <- as.factor(rep(seq(300,600,50),2)); DepCat
facttonum(DepCat)

## Writing Functions
# R-chunk 6  page 25
#Outline of a function's structure

functionname <- function(argument1, fun,...) {
  # body of the function
  #
  # the input arguments and body of a function can include other
  # functions, which may have their own arguments, which is what
  # the ... is for. One can include other inputs that are used but
  # not defined early on and may depend on what function is brought
  # into the main function. See for example negLL(), and others
  answer <- fun(argument1) + 2
  return(answer)
} # end of functionname
functionname(c(1,2,3,4,5),mean)  # = mean(1,2,3,4,5)= 3 + 2 = 5

### Simple Functions
# R-chunk 7  page 26
# Implement the von Bertalanffy curve in multiple ways

ages <- 1:20
nages <- length(ages)
Linf <- 50;  K <- 0.2;  t0 <- -0.75
# first try a for loop to calculate length for each age
loopLt <- numeric(nages)
for (ag in ages) loopLt[ag] <- Linf * (1 - exp(-K * (ag - t0)))
# the equations are automatically vectorized so more efficient
vecLt <- Linf * (1 - exp(-K * (ages - t0))) # or we can convert
# the equation into a function and use it again and again
vB <- function(pars,inages) { # requires pars=c(Linf,K,t0)
  Lt <- pars[1] * (1 - exp(-pars[2] * (inages - pars[3])))
  return(Lt)
}
funLt <- vB(c(Linf,K,t0),ages)
ans <- cbind(ages,funLt,vecLt,loopLt)


# R-chunk 8  page 26 - code not shown in book.
# Tabulate the ans from chunk 7

library(knitr)  # needed for the function knitr - pretty tables
kable(halftable(ans,yearcol="ages",subdiv=2),digits=c(0,3,3,3,0,3,3,3))


# R-chunk 9  page 27
#A vB function with some input error checking

vB <- function(pars,inages) { # requires pars=c(Linf,K,t0)
  if (is.numeric(pars) & is.numeric(inages)) {
    Lt <- pars[1] * (1 - exp(-pars[2] * (inages - pars[3])))
  } else { stop(cat("Not all input values are numeric! \n")) }
  return(Lt)
}
param <- c(50, 0.2,"-0.75")
funLt <- vB(as.numeric(param),ages) #try without the as.numeric
halftable(cbind(ages,funLt))


### Scoping of Objects
# R-chunk 10   page 29
# demonstration that the globel environment is 'visible' inside a
# a function it calls, but the function's environment remains
# invisible to the global or calling environment

vBscope <- function(pars) { # requires pars=c(Linf,K,t0)
  rhside <- (1 - exp(-pars[2] * (ages - pars[3])))
  Lt <- pars[1] * rhside
  return(Lt)
}
ages <- 1:10; param <- c(50,0.2,-0.75)
vBscope(param)
try(rhside)    # note the use of try() which can trap errors ?try


### Function Inputs and Outputs
# R-chunk 11  page 30
#Bring the data-set schaef into the working of global environment

data(schaef)


# R-chunk 12   page 30  Table 2.2 code not shown
#Tabulate the data held in schaef. Needs knitr

kable(halftable(schaef,yearcol="year",subdiv=2),digits=c(0,0,0,4))

# R-chunk 13  page 30
#examine the properties of the data-set schaef

class(schaef)
a <- schaef[1:5,2]
b <- schaef[1:5,"catch"]
c <- schaef$catch[1:5]
cbind(a,b,c)
mschaef <- as.matrix(schaef)
mschaef[1:5,"catch"]  # ok
d <- try(mschaef$catch[1:5]) #invalid for matrices
d  # had we not used try()eveerything would have stopped.


# R-chunk 14  page 31
#Convert column names of a data.frame or matrix to lowercase

dolittle <- function(indat) {
  indat1 <- as.data.frame(indat)
  colnames(indat) <- tolower(colnames(indat))
  return(list(dfdata=indat1,indat=as.matrix(indat)))
} # return the original and the new version
colnames(schaef) <- toupper(colnames(schaef))
out <- dolittle(schaef)
str(out, width=63, strict.width="cut")


# R-chunk 15  page 32
#Could have used an S3 plot method had we defined a class   Fig.2.2

plotspmdat(schaef) # examine the code as an eg of a custom plot

## End(Not run)

Description

chapter3 is not an active function but rather acts as a repository for the various example code chunks found in chapter3. There are 27 r-code chunks in chapter3.

Usage

chapter3(verbose = TRUE)

Arguments

Examples

## Not run: 
### The Discrete Logistic Model    
# R-chunk 1 page 36
# Code to produce Figure 3.1. Note the two one-line functions 
   
surprod <- function(Nt,r,K) return((r*Nt)*(1-(Nt/K)))    
densdep <- function(Nt,K) return((1-(Nt/K)))    
r <- 1.2; K <- 1000.0; Nt <- seq(10,1000,10)  
oldpar <- par(no.readonly=TRUE) # this line not in book
# plotprep(width=7, height=5, newdev=FALSE)  
par(mfrow=c(2,1),mai=c(0.4,0.4,0.05,0.05),oma=c(0.0,0,0.0,0.0))     
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)     
plot1(Nt,surprod(Nt,r,K),xlab="Population Nt",defpar=FALSE,    
      ylab="Production")    
plot1(Nt,densdep(Nt,K),xlab="Population Nt",defpar=FALSE,    
      ylab="Density-Dependence")    
par(oldpar)  # this line not in book


### Dynamic Behaviour    
# R-chunk 2 page 38
#Code for Figure 3.2. Try varying the value of rv from 0.5-2.8 
   
yrs <- 100; rv=2.8;  Kv <- 1000.0; Nz=100; catch=0.0; p=1.0    
ans <- discretelogistic(r=rv,K=Kv,N0=Nz,Ct=catch,Yrs=yrs,p=p)    
avcatch <- mean(ans[(yrs-50):yrs,"nt"],na.rm=TRUE) #used in text    
label <- paste0("r=",rv," K=",Kv," Ct=",catch, " N0=",Nz," p=",p=p)   
oldpar <- par(no.readonly=TRUE) # this line not in book
plot(ans, main=label, cex=0.9, font=7) #Schaefer dynamics    
par(oldpar)   # this line not in book


# R-chunk 3 page 39
#run discrete logistic dynamics for 600 years  
  
yrs=600    
ans <- discretelogistic(r=2.55,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)    

# R-chunk 4  page 40, code not in the book
#tabulate the last 30 years of the dynamics   needs knitr
library(knitr) 
kable(halftable(ans[(yrs-29):yrs,],yearcol="year",subdiv=3),digits=c(0,1,1,0,1,1,0,1,1))    


### Finding Boundaries between Behaviours.    
# R-chunk 5 page 40
#run discretelogistic and search for repeated values of Nt    

yrs <- 600    
ans <- discretelogistic(r=2.55,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)    
avt <- round(apply(ans[(yrs-100):(yrs-1),2:3],1,mean),2)    
count <- table(avt)    
count[count > 1] # with r=2.55 you should find an 8-cycle limit    

# R-chunk 6  page 41
#searches for unique solutions given an r value  see Table 3.2  

testseq <- seq(1.9,2.59,0.01)    
nseq <- length(testseq)    
result <- matrix(0,nrow=nseq,ncol=2,    
                 dimnames=list(testseq,c("r","Unique")))    
yrs <- 600    
for (i in 1:nseq) {  # i = 31    
  rval <- testseq[i]    
  ans <- discretelogistic(r=rval,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)    
  ans <- ans[-yrs,] # remove last year, see str(ans) for why    
  ans[,"nt1"] <- round(ans[,"nt1"],3) #try hashing this out    
  result[i,] <- c(rval,length(unique(tail(ans[,"nt1"],100))))    
}    


# R-chunk 7  page 41 - 42, Table 3.2. Code not in the book.
#unique repeated Nt values 100 = non-equilibrium or chaos   
 
kable(halftable(result,yearcol = "r"),)    


### Classical Bifurcation Diagram of Chaos    
# R-chunk 8 pages 42 - 43
#the R code for the bifurcation function   
 
bifurcation <- function(testseq,taill=100,yrs=1000,limy=0,incx=0.001){    
  nseq <- length(testseq)    
  result <- matrix(0,nrow=nseq,ncol=2,    
                   dimnames=list(testseq,c("r","Unique Values")))    
  result2 <- matrix(NA,nrow=nseq,ncol=taill)    
  for (i in 1:nseq) {      
    rval <- testseq[i]    
    ans <- discretelogistic(r=rval,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)    
    ans[,"nt1"] <- round(ans[,"nt1"],4)    
    result[i,] <- c(rval,length(unique(tail(ans[,"nt1"],taill))))    
    result2[i,] <- tail(ans[,"nt1"],taill)    
  }      
  if (limy[1] == 0) limy <- c(0,getmax(result2,mult=1.02))  
    
  oldpar <- parset() #plot taill values against taill of each r value  
  on.exit(par(oldpar))    #  this line not in book 
  plot(rep(testseq[1],taill),result2[1,],type="p",pch=16,cex=0.1,    
       ylim=limy,xlim=c(min(testseq)*(1-incx),max(testseq)*(1+incx)),    
       xlab="r value",yaxs="i",xaxs="i",ylab="Equilibrium Numbers",    
       panel.first=grid())    
  for (i in 2:nseq)    
    points(rep(testseq[i],taill),result2[i,],pch=16,cex=0.1)    
  return(invisible(list(result=result,result2=result2)))    
} # end of bifurcation    


# R-chunk 9 page 43 
#Alternative r value arrangements for you to try; Fig 3.3    
#testseq <- seq(2.847,2.855,0.00001) #hash/unhash as needed    
#bifurcation(testseq,limy=c(600,740),incx=0.0001) # t    
#testseq <- seq(2.6225,2.6375,0.00001) # then explore     
#bifurcation(testseq,limy=c(660,730),incx=0.0001)  
   
testseq <- seq(1.9,2.975,0.0005) # modify to explore    
bifurcation(testseq,limy=0)      


### The Effect of Fishing on Dynamics    
# R-chunk 10  page 43 - 44.
#Effect of catches on stability properties of discretelogistic   
 
yrs=50; Kval=1000.0    
nocatch <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=0,Yrs=yrs)    
catch50 <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=50,Yrs=yrs)    
catch200 <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=200,Yrs=yrs)    
catch300 <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=300,Yrs=yrs)    


# R-chunk 11 page 45
#Effect of different catches on n-cyclic behaviour Fig3.4   
 
plottime <- function(x,ylab) {    
  yrs <- nrow(x)    
  plot1(x[,"year"],x[,"nt"],ylab=ylab,defpar=FALSE)    
  avB <- round(mean(x[(yrs-40):yrs,"nt"],na.rm=TRUE),3)    
  mtext(avB,side=1,outer=F,line=-1.1,font=7,cex=1.0)     
} # end of plottime    
#the oma argument is used to adjust the space around the graph 
oldpar <- par(no.readonly=TRUE) # this line not in book   
par(mfrow=c(2,2),mai=c(0.25,0.4,0.05,0.05),oma=c(1.0,0,0.25,0))     
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)      
plottime(nocatch,"Catch = 0")    
plottime(catch50,"Catch = 50")    
plottime(catch200,"Catch = 200")    
plottime(catch300,"Catch = 300")    
mtext("years",side=1,outer=TRUE,line=-0.2,font=7,cex=1.0)     
par(oldpar)

# R-chunk 12 page 46
#Phase plot for Schaefer model Fig 3.5    

plotphase <- function(x,label,ymax=0) { #x from discretelogistic    
  yrs <- nrow(x)    
  colnames(x) <- tolower(colnames(x))    
  if (ymax[1] == 0) ymax <- getmax(x[,c(2:3)])    
  plot(x[,"nt"],x[,"nt1"],type="p",pch=16,cex=1.0,ylim=c(0,ymax),    
       yaxs="i",xlim=c(0,ymax),xaxs="i",ylab="nt1",xlab="",    
       panel.first=grid(),col="darkgrey")    
  begin <- trunc(yrs * 0.6) #last 40% of yrs = 20, when yrs=50    
  points(x[begin:yrs,"nt"],x[begin:yrs,"nt1"],pch=18,col=1,cex=1.2)    
  mtext(label,side=1,outer=F,line=-1.1,font=7,cex=1.2)     
} # end of plotphase    
oldpar <- par(no.readonly=TRUE) # this line not in book   
par(mfrow=c(2,2),mai=c(0.25,0.25,0.05,0.05),oma=c(1.0,1.0,0,0))     
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)      
plotphase(nocatch,"Catch = 0",ymax=1300)    
plotphase(catch50,"Catch = 50",ymax=1300)    
plotphase(catch200,"Catch = 200",ymax=1300)    
plotphase(catch300,"Catch = 300",ymax=1300)    
mtext("nt",side=1,outer=T,line=0.0,font=7,cex=1.0)    
mtext("nt+1",side=2,outer=T,line=0.0,font=7,cex=1.0)    
par(oldpar)    # this line not in book

### Determinism    
## Age-Structured Modelling Concepts    
### Survivorship in a Cohort    
# R-chunk 13 pages 48 - 49
#Exponential population declines under different Z. Fig 3.6   
 
yrs <- 50;  yrs1 <- yrs + 1 # to leave room for B[0]    
years <- seq(0,yrs,1)    
B0 <- 1000        # now alternative total mortality rates    
Z <- c(0.05,0.1,0.2,0.4,0.55)     
nZ <- length(Z)    
Bt <- matrix(0,nrow=yrs1,ncol=nZ,dimnames=list(years,Z))    
Bt[1,] <- B0    
for (j in 1:nZ) for (i in 2:yrs1) Bt[i,j] <- Bt[(i-1),j]*exp(-Z[j])    
oldp <- plot1(years,Bt[,1],xlab="Years",ylab="Population Size",lwd=2)    
if (nZ > 1) for (j in 2:nZ) lines(years,Bt[,j],lwd=2,col=j,lty=j)    
legend("topright",legend=paste0("Z = ",Z),col=1:nZ,lwd=3,    
       bty="n",cex=1,lty=1:5)     
par(oldp)  # this line not in book

### Instantaneous vs Annual Mortality Rates    
# R-chunk 14 page 51
#Prepare matrix of harvest rate vs time to appoximate F   
 
Z <- -log(0.5)    
timediv <- c(2,4,12,52,365,730,2920,8760,525600)    
yrfrac <- 1/timediv    
names(yrfrac) <- c("6mth","3mth","1mth","1wk","1d","12h","3h","1h","1m")    
nfrac <- length(yrfrac)    
columns <- c("yrfrac","divisor","yrfracH","Remain")    
result <- matrix(0,nrow=nfrac,ncol=length(columns),    
                 dimnames=list(names(yrfrac),columns))    
for (i in 1:nfrac) {    
  timestepmort <- Z/timediv[i]     
  N <- 1000    
  for (j in 1:timediv[i]) N <- N * (1-timestepmort)    
  result[i,] <- c(yrfrac[i],timediv[i],timestepmort,N)    
}    


# R-chunk 15 page 51 Table 3.3, code not shown in book
#output of constant Z for shorter and shorter periods    

kable(result,digits=c(10,0,8,4))    


# R-chunk 16 page 51
#Annual harvest rate against instantaneous F, Fig 3.7  
  
Fi <- seq(0.001,2,0.001)    
H <- 1 - exp(-Fi)    
oldpar <- parset()  # a wrapper for simplifying defining the par values    
plot(Fi,H,type="l",lwd=2,panel.first=grid(),xlab="Instantaneous Fishing Mortality F",    
     ylab="Annual Proportion Mortality H")    
lines(c(0,1),c(0,1),lwd=2,lty=2,col=2)    
par(oldpar)   # this line not in book


## Simple Yield per Recruit    
# R-chunk 17  page 53
# Simple Yield-per-Recruit see Russell (1942)   
 
age <- 1:11;  nage <- length(age); N0 <- 1000  # some definitions    
# weight-at-age values    
WaA <- c(NA,0.082,0.175,0.283,0.4,0.523,0.7,0.85,0.925,0.99,1.0)    
# now the harvest rates    
H <- c(0.01,0.06,0.11,0.16,0.21,0.26,0.31,0.36,0.55,0.8)    
nH <- length(H)    
NaA <- matrix(0,nrow=nage,ncol=nH,dimnames=list(age,H)) # storage    
CatchN <- NaA;  CatchW <- NaA      # define some storage matrices    
for (i in 1:nH) {                # loop through the harvest rates    
  NaA[1,i] <- N0  # start each harvest rate with initial numbers    
  for (age in 2:nage) {  # loop through over-simplified dynamics    
    NaA[age,i] <- NaA[(age-1),i] * (1 - H[i])    
    CatchN[age,i] <- NaA[(age-1),i] - NaA[age,i]    
  }    
  CatchW[,i] <- CatchN[,i] * WaA    
}                      # transpose the vector of total catches to    
totC <- t(colSums(CatchW,na.rm=TRUE))   # simplify later printing    


# R-chunk 18 page 54 Table 3.4 code not shown in book
#Tabulate numbers-at-age for different harvest rates  needs knitr
  
kable(NaA,digits=c(0,0,0,0,0,0,0,0,1,1),row.names=TRUE)    


# R-chunk 19 page 54, Table 3.5, code not shown in book.
#Tabulate Weight-at-age for different harvest rates   
 
kable(CatchW[2:11,],digits=c(2,2,2,2,2,2,2,2,2,2),row.names=TRUE)    


# R-chunk 20 page 54, Table 3.6, code not shown in book.
#Total weights vs Harvest rate   
 
kable(totC,digits=c(1,1,1,1,1,1,1,1,1,1))    


# R-chunk 21 page 55
#Use MQMF::plot1 for a quick plot of the total catches. Figure 3.8    

oldpar <- plot1(H,totC,xlab="Harvest Rate",ylab="Total Yield",lwd=2)    
par(oldpar) # to reset the par values if desired

### Selectivity in Yield-per-Recruit    
# R-chunk 22 Page 56
#Logistic S shaped cureve for maturity    

ages <- seq(0,50,1)    
sel1 <- mature(-3.650425,0.146017,sizeage=ages) #-3.65/0.146=25    
sel2 <- mature(-6,0.2,ages)    
sel3 <- mature(-6,0.24,ages)    
oldp <- plot1(ages,sel1,xlab="Age Yrs",ylab="Selectivity",cex=0.75,lwd=2)    
lines(ages,sel2,col=2,lwd=2,lty=2)    
lines(ages,sel3,col=3,lwd=2,lty=3)    
abline(v=25,col="grey",lty=2)     
abline(h=c(0.25,0.5,0.75),col="grey",lty=2)    
legend("topleft",c("25_15.04","30_10.986","25_9.155"),col=c(1,2,3),    
       lwd=3,cex=1.1,bty="n",lty=1:3)    
par(oldp)

### The Baranov Catch Equation    
# R-chunk 23 Page 58
# Baranov catch equation  
  
age <- 0:12;  nage <- length(age)     
sa <-mature(-4,2,age) #selectivity-at-age    
H <- 0.2;  M <- 0.35    
FF <- -log(1 - H)#Fully selected instantaneous fishing mortality    
Ft <- sa * FF     # instantaneous Fishing mortality-at-age    
N0 <- 1000    
out <- cbind(bce(M,Ft,N0,age),"Select"=sa)  # out becomes Table 3.7    


# R-chunk 24 page 59, Table 3.7, code not shown in book.
#tabulate output from Baranov Catch Equations     

kable(out,digits=c(3,3,3,3))    


### Growth and Weight-at-Age   
## Full Yield-per-Recruit    
# R-chunk 25 Page 60 - 61
# A more complete YPR analysis    

age <- 0:20;  nage <- length(age) #storage vectors and matrices    
laa <- vB(c(50.0,0.25,-1.5),age) # length-at-age    
WaA <- (0.015 * laa ^ 3.0)/1000  # weight-at-age as kg    
H <- seq(0.01,0.65,0.05);  nH <- length(H)       
FF <- round(-log(1 - H),5)  # Fully selected fishing mortality    
N0 <- 1000    
M <- 0.1    
numt <- matrix(0,nrow=nage,ncol=nH,dimnames=list(age,FF))    
catchN <- matrix(0,nrow=nage,ncol=nH,dimnames=list(age,FF))    
as50 <- c(1,2,3)      
yield <- matrix(0,nrow=nH,ncol=length(as50),dimnames=list(H,as50))    
for (sel in 1:length(as50)) {    
  sa <- logist(as50[sel],1.0,age)  # selectivity-at-age    
  for (harv in 1:nH) {    
    Ft <- sa * FF[harv]      # Fishing mortality-at-age    
    out <- bce(M,Ft,N0,age)    
    numt[,harv] <- out[,"Nt"]    
    catchN[,harv] <- out[,"Catch"]    
    yield[harv,sel] <- sum(out[,"Catch"] * WaA,na.rm=TRUE)    
  } # end of harv loop    
} # end of sel loop    


# R-chunk 26  Page 61
#A full YPR analysis  Figure 3.10    

oldp <- plot1(H,yield[,3],xlab="Harvest Rate",ylab="Yield",cex=0.75,lwd=2)    
lines(H,yield[,2],lwd=2,col=2,lty=2)    
lines(H,yield[,1],lwd=2,col=3,lty=3)    
legend("bottomright",legend=as50,col=c(3,2,1),lwd=3,bty="n",    
       cex=1.0,lty=c(3,2,1))     
par(oldp)

# R-chunk 27 page 62, Table 3.8, code not shown in book.
#Tabulate yield-per-recruit using Baranoc catch equation   
 
kable(yield,digits=c(2,3,3,3))    


## End(Not run)

Description

chapter4 is not an active function but rather acts as a repository for the various example code chunks found in chapter4. There are 47 r-code chunks in chapter3.

Usage

chapter4(verbose = TRUE)

Arguments

Examples

## Not run: 
# All the example code from  # Model Parameter Estimation    
# Model Parameter Estimation    
## Introduction    
### Optimization    
## Criteria of Best Fit    
## Model Fitting in R    
### Model Requirements    
### A Length-at-Age Example    
### Alternative Models of Growth    
## Sum of Squared Residual Deviations    
### Assumptions of Least-Squares    
### Numerical Solutions   
 
# R-chunk 1  Page 75 
#setup optimization using growth and ssq    
#convert equations 4.4 to 4.6 into vectorized R functions    
#These will over-write the same functions in the MQMF package    

data(LatA)      # try ?LatA   assumes library(MQMF) already run    
vB <- function(p, ages) return(p[1]*(1-exp(-p[2]*(ages-p[3]))))    
Gz <- function(p, ages) return(p[1]*exp(-p[2]*exp(p[3]*ages)))    
mm <- function(p, ages) return((p[1]*ages)/(p[2] + ages^p[3]))    
#specific function to calc ssq. The ssq within MQMF is more    
ssq <- function(p,funk,agedata,observed) {        #general and is    
  predval <- funk(p,agedata)        #not limited to p and agedata    
  return(sum((observed - predval)^2,na.rm=TRUE))    
} #end of ssq     
# guess starting values for Linf, K, and t0, names not needed    
pars <- c("Linf"=27.0,"K"=0.15,"t0"=-2.0) #ssq should=1478.449    
ssq(p=pars, funk=vB, agedata=LatA$age, observed=LatA$length)     
# try misspelling LatA$Length with a capital. What happens?    

### Passing Functions as Arguments to other Functions    
# R-chunk 2  Page 76 
# Illustrates use of names within function arguments    

vB <- function(p,ages) return(p[1]*(1-exp(-p[2] *(ages-p[3]))))    
ssq <- function(funk,observed,...) { # only define ssq arguments    
  predval <- funk(...) # funks arguments are implicit    
  return(sum((observed - predval)^2,na.rm=TRUE))    
} # end of ssq     
pars <- c("Linf"=27.0,"K"=0.15,"t0"=-2.0) # ssq should = 1478.449    
ssq(p=pars, funk=vB, ages=LatA$age, observed=LatA$length) #if no    
ssq(vB,LatA$length,pars,LatA$age) # name order is now vital!    


# R-chunk 3  Page 77
# Illustrate a problem with calling a function in a function    
# LatA$age is typed as LatA$Age but no error, and result = 0    

ssq(funk=vB, observed=LatA$length, p=pars, ages=LatA$Age) # !!!    

### Fitting the Models    
# R-chunk 4  Page 77
#plot the LatA data set   Figure 4.2    

oldpar <- parset()   # parset and getmax are two MQMF functions     
ymax <- getmax(LatA$length) # simplifies use of base graphics. For    
# full colour, with the rgb as set-up below, there must be >= 5 obs    
plot(LatA$age,LatA$length,type="p",pch=16,cex=1.2,xlab="Age Years",     
     ylab="Length cm",col=rgb(1,0,0,1/5),ylim=c(0,ymax),yaxs="i",    
     xlim=c(0,44),panel.first=grid()) 
par(oldpar) # this line not in book   

# R-chunk 5  Pages 78 - 79 
# use nlm to fit 3 growth curves to LatA, only p and funk change   

ages <- 1:max(LatA$age) # used in comparisons     
pars <- c(27.0,0.15,-2.0) # von Bertalanffy    
bestvB <- nlm(f=ssq,funk=vB,observed=LatA$length,p=pars,    
              ages=LatA$age,typsize=magnitude(pars))    
outfit(bestvB,backtran=FALSE,title="vB"); cat("\n")     
pars <- c(26.0,0.7,-0.5) # Gompertz    
bestGz <- nlm(f=ssq,funk=Gz,observed=LatA$length,p=pars,    
              ages=LatA$age,typsize=magnitude(pars))    
outfit(bestGz,backtran=FALSE,title="Gz"); cat("\n")     
pars <- c(26.2,1.0,1.0) # Michaelis-Menton - first start point    
bestMM1 <- nlm(f=ssq,funk=mm,observed=LatA$length,p=pars,    
               ages=LatA$age,typsize=magnitude(pars))    
outfit(bestMM1,backtran=FALSE,title="MM"); cat("\n")    
pars <- c(23.0,1.0,1.0) # Michaelis-Menton - second start point    
bestMM2 <- nlm(f=ssq,funk=mm,observed=LatA$length,p=pars,    
               ages=LatA$age,typsize=magnitude(pars))    
outfit(bestMM2,backtran=FALSE,title="MM2"); cat("\n")     

# R-chunk 6  Page 81 
#The use of args() and formals()     

args(nlm) # formals(nlm) uses more screen space. Try yourself.    


# R-chunk 7  Page 81, code not in the book 
#replacement for args(nlm) to keep within page borders without truncation   

{cat("function (f, p, ..., hessian = FALSE, typsize = rep(1,\n")    
  cat("  length(p)),fscale = 1, print.level = 0, ndigit = 12, \n")    
  cat("  gradtol = 1e-06, stepmax = max(1000 * \n")    
  cat("  sqrt(sum((p/typsize)^2)), 1000), steptol = 1e-06, \n")     
  cat("  iterlim = 100, check.analyticals = TRUE)\n")}    


# R-chunk 8  Pages 81 - 82 
#Female length-at-age + 3 growth fitted curves Figure 4.3    

predvB <- vB(bestvB$estimate,ages) #get optimumpredicted lengths    
predGz <- Gz(bestGz$estimate,ages) # using the outputs    
predmm <- mm(bestMM2$estimate,ages) #from the nlm analysis above    
ymax <- getmax(LatA$length) #try ?getmax or getmax [no brackets]    
xmax <- getmax(LatA$age)  #there is also a getmin, not used here    
oldpar <- parset(font=7) #or use parsyn() to prompt for par syntax    
plot(LatA$age,LatA$length,type="p",pch=16, col=rgb(1,0,0,1/5),    
     cex=1.2,xlim=c(0,xmax),ylim=c(0,ymax),yaxs="i",xlab="Age",    
     ylab="Length (cm)",panel.first=grid())    
lines(ages,predvB,lwd=2,col=4)        # vB    col=4=blue    
lines(ages,predGz,lwd=2,col=1,lty=2)  # Gompertz  1=black    
lines(ages,predmm,lwd=2,col=3,lty=3)  # MM        3=green    
#notice the legend function and its syntax.    
legend("bottomright",cex=1.2,c("von Bertalanffy","Gompertz",    
                   "Michaelis-Menton"),col=c(4,1,3),lty=c(1,2,3),lwd=3,bty="n") 
par(oldpar)  # this line not in book                       

### Objective Model Selection    
### The Influence of Residual Error Choice on Model Fit    
# R-chunk 9  Page 84 - 85 
# von Bertalanffy     

pars <- c(27.25,0.15,-3.0)    
bestvBN <- nlm(f=ssq,funk=vB,observed=LatA$length,p=pars,    
               ages=LatA$age,typsize=magnitude(pars),iterlim=1000)    
outfit(bestvBN,backtran=FALSE,title="Normal errors"); cat("\n")     
# modify ssq to account for log-normal errors in ssqL    
ssqL <- function(funk,observed,...) {    
  predval <- funk(...)    
  return(sum((log(observed) - log(predval))^2,na.rm=TRUE))    
} # end of ssqL    
bestvBLN <- nlm(f=ssqL,funk=vB,observed=LatA$length,p=pars,    
                ages=LatA$age,typsize=magnitude(pars),iterlim=1000)    
outfit(bestvBLN,backtran=FALSE,title="Log-Normal errors")    


# R-chunk 10  Pages 85 - 86
# Now plot the resultibng two curves and the data Fig 4.4    

predvBN <- vB(bestvBN$estimate,ages)     
predvBLN <- vB(bestvBLN$estimate,ages)     
ymax <- getmax(LatA$length)     
xmax <- getmax(LatA$age)        
oldpar <- parset()                  
plot(LatA$age,LatA$length,type="p",pch=16, col=rgb(1,0,0,1/5),    
     cex=1.2,xlim=c(0,xmax),ylim=c(0,ymax),yaxs="i",xlab="Age",    
     ylab="Length (cm)",panel.first=grid())    
lines(ages,predvBN,lwd=2,col=4,lty=2)   # add Normal dashed    
lines(ages,predvBLN,lwd=2,col=1)        # add Log-Normal solid    
legend("bottomright",c("Normal Errors","Log-Normal Errors"),    
       col=c(4,1),lty=c(2,1),lwd=3,bty="n",cex=1.2)    
par(oldpar)

### Remarks on Initial Model Fitting    
## Maximum Likelihood     
### Introductory Examples    
# R-chunk 11  Page 88
# Illustrate Normal random likelihoods. see Table 4.1    

set.seed(12345)       # make the use of random numbers repeatable    
x <- rnorm(10,mean=5.0,sd=1.0)      # pseudo-randomly generate 10     
avx <- mean(x)                      # normally distributed values    
sdx <- sd(x)          # estimate the mean and stdev of the sample              
L1 <- dnorm(x,mean=5.0,sd=1.0)   # obtain likelihoods, L1, L2 for     
L2 <- dnorm(x,mean=avx,sd=sdx)    # each data point for both sets    
result <- cbind(x,L1,L2,"L2gtL1"=(L2>L1))      # which is larger?    
result <- rbind(result,c(NA,prod(L1),prod(L2),1)) # result+totals    
rownames(result) <- c(1:10,"product")    
colnames(result) <- c("x","original","estimated","est > orig")    


# R-chunk 12  page 88, Table 4.1, code not shown in book.
#tabulate results of Normal Likelihoods    

kable(result,digits=c(4,8,8,0),row.names = TRUE, caption='(ref:tab401)')    

# R-chunk 13  Page 89 
# some examples of pnorm, dnorm, and qnorm, all mean = 0    

cat("x = 0.0        Likelihood =",dnorm(0.0,mean=0,sd=1),"\n")     
cat("x = 1.95996395 Likelihood =",dnorm(1.95996395,mean=0,sd=1),"\n")     
cat("x =-1.95996395 Likelihood =",dnorm(-1.95996395,mean=0,sd=1),"\n")     
# 0.5 = half cumulative distribution    
cat("x = 0.0        cdf = ",pnorm(0,mean=0,sd=1),"\n")     
cat("x = 0.6744899  cdf = ",pnorm(0.6744899,mean=0,sd=1),"\n")    
cat("x = 0.75       Quantile =",qnorm(0.75),"\n") # reverse pnorm    
cat("x = 1.95996395 cdf = ",pnorm(1.95996395,mean=0,sd=1),"\n")    
cat("x =-1.95996395 cdf = ",pnorm(-1.95996395,mean=0,sd=1),"\n")    
cat("x = 0.975      Quantile =",qnorm(0.975),"\n") # expect ~1.96    
# try x <- seq(-5,5,0.2); round(dnorm(x,mean=0.0,sd=1.0),5)    

## Likelihoods from the Normal Distribution    
# R-chunk 14   Page 90
# Density plot and cumulative distribution for Normal   Fig 4.5    

x <- seq(-5,5,0.1)  # a sequence of values around a mean of 0.0    
NL <- dnorm(x,mean=0,sd=1.0)   # normal likelihoods for each X    
CD <- pnorm(x,mean=0,sd=1.0)   # cumulative density vs X    
oldp <- plot1(x,CD,xlab="x = StDev from Mean",ylab="Likelihood and CDF")    
lines(x,NL,lwd=3,col=2,lty=3) # dashed line as these are points    
abline(h=0.5,col=4,lwd=1)  
par(oldp)  

# R-chunk 15  Pages 91 - 92 
#function facilitates exploring different polygons Fig 4.6    

plotpoly <- function(mid,delta,av=5.0,stdev=1.0) {    
  neg <- mid-delta;  pos <- mid+delta    
  pdval <- dnorm(c(mid,neg,pos),mean=av,sd=stdev)    
  polygon(c(neg,neg,mid,neg),c(pdval[2],pdval[1],pdval[1],    
                               pdval[2]),col=rgb(0.25,0.25,0.25,0.5))    
  polygon(c(pos,pos,mid,pos),c(pdval[1],pdval[3],pdval[1],    
                               pdval[1]),col=rgb(0,1,0,0.5))       
  polygon(c(mid,neg,neg,mid,mid),    
          c(0,0,pdval[1],pdval[1],0),lwd=2,lty=1,border=2)    
  polygon(c(mid,pos,pos,mid,mid),    
          c(0,0,pdval[1],pdval[1],0),lwd=2,lty=1,border=2)     
  text(3.395,0.025,paste0("~",round((2*(delta*pdval[1])),7)),  
       cex=1.1,pos=4)    
  return(2*(delta*pdval[1])) # approx probability, see below    
} # end of plotpoly, a temporary function to enable flexibility    
#This code can be re-run with different values for delta    
x <- seq(3.4,3.6,0.05) # where under the normal curve to examine    
pd <- dnorm(x,mean=5.0,sd=1.0) #prob density for each X value    
mid <- mean(x)        
delta <- 0.05  # how wide either side of the sample mean to go?     
oldpar <- parset() #pre-defined MQMF base graphics set-up for par    
ymax <- getmax(pd) # find maximum y value for the plot    
plot(x,pd,type="l",xlab="Variable x",ylab="Probability Density",    
     ylim=c(0,ymax),yaxs="i",lwd=2,panel.first=grid())    
approxprob <- plotpoly(mid,delta)  #use function defined above 
par(oldpar)  # this line not in the book

### Equivalence with Sum-of-Squares     
### Fitting a Model to Data using Normal Likelihoods    
# R-chunk 16  Page 94 
#plot of length-at-age data  Fig 4.7    

data(LatA) # load the redfish data set into memory and plot it    
ages <- LatA$age;  lengths <- LatA$length    
oldpar <- plot1(ages,lengths,xlab="Age",ylab="Length",type="p",cex=0.8,    
      pch=16,col=rgb(1,0,0,1/5))    
par(oldpar)


# R-chunk 17  Page 95
# Fit the vB growth curve using maximum likelihood    

pars <- c(Linf=27.0,K=0.15,t0=-3.0,sigma=2.5) # starting values    
# note, estimate for sigma is required for maximum likelihood    
ansvB <- nlm(f=negNLL,p=pars,funk=vB,observed=lengths,ages=ages,    
             typsize=magnitude(pars))    
outfit(ansvB,backtran=FALSE,title="vB by minimum -veLL")    

# R-chunk 18 Page 96
#Now fit the Michaelis-Menton curve    

pars <- c(a=23.0,b=1.0,c=1.0,sigma=3.0) # Michaelis-Menton  
ansMM <- nlm(f=negNLL,p=pars,funk=mm,observed=lengths,ages=ages,    
             typsize=magnitude(pars))    
outfit(ansMM,backtran=FALSE,title="MM by minimum -veLL")    

# R-chunk 19  Page 96 
#plot optimum solutions for vB and mm. Fig 4.8    

Age <- 1:max(ages) # used in comparisons     
predvB <- vB(ansvB$estimate,Age) #optimum solution    
predMM <- mm(ansMM$estimate,Age) #optimum solution    
oldpar <- parset()               # plot the deata points first  
plot(ages,lengths,xlab="Age",ylab="Length",type="p",pch=16,    
     ylim=c(10,33),panel.first=grid(),col=rgb(1,0,0,1/3))    
lines(Age,predvB,lwd=2,col=4)     # then add the growth curves  
lines(Age,predMM,lwd=2,col=1,lty=2)    
legend("bottomright",c("von Bertalanffy","Michaelis-Menton"),    
       col=c(4,1),lwd=3,bty="n",cex=1.2,lty=c(1,2)) 
par(oldpar)   

# R-chunk 20   Pages 96 - 97
# residual plot for vB curve   Fig 4.9    

predvB <- vB(ansvB$estimate,ages) # predicted values for age data    
resids <- lengths - predvB               # calculate vB residuals     
oldpar <- plot1(ages,resids,type="p",col=rgb(1,0,0,1/3),
           xlim=c(0,43),pch=16,xlab="Ages Years",ylab="Residuals")    
abline(h=0.0,col=1,lty=2)    # emphasize the zero line  
par(oldpar)

## Log-Normal Likelihoods     
### Simplification of Log-Normal Likelihoods    
### Log-Normal Properties    
# R-chunk 21  Page 100 
# meanlog and sdlog affects on mode and spread of lognormal Fig 4.10     

x <- seq(0.05,5.0,0.01)  # values must be greater than 0.0    
y <- dlnorm(x,meanlog=0,sdlog=1.2,log=FALSE) #dlnorm=likelihoods    
y2 <- dlnorm(x,meanlog=0,sdlog=1.0,log=FALSE)#from log-normal     
y3 <- dlnorm(x,meanlog=0,sdlog=0.6,log=FALSE)#distribution     
y4 <- dlnorm(x,0.75,0.6)         #log=TRUE = log-likelihoods    
oldpar <- parset(plots=c(1,2)) #MQMF base plot formatting function    
plot(x,y3,type="l",lwd=2,panel.first=grid(),    
     ylab="Log-Normal Likelihood")    
lines(x,y,lwd=2,col=2,lty=2)    
lines(x,y2,lwd=2,col=3,lty=3)    
lines(x,y4,lwd=2,col=4,lty=4)    
legend("topright",c("meanlog sdlog","    0.0      0.6    0.0",    
                    "      1.0","    0.0      1.2","    0.75    0.6"),    
       col=c(0,1,3,2,4),lwd=3,bty="n",cex=1.0,lty=c(0,1,3,2,4))    
plot(log(x),y3,type="l",lwd=2,panel.first=grid(),ylab="")    
lines(log(x),y,lwd=2,col=2,lty=2)    
lines(log(x),y2,lwd=2,col=3,lty=3)    
lines(log(x),y4,lwd=2,col=4,lty=4)  
par(oldpar)  # return par to old settings; this line not in book


# R-chunk 22  Pages 100 - 101

set.seed(12354) # plot random log-normal numbers as Fig 4.11    
meanL <- 0.7;   sdL <- 0.5  # generate 5000 random log-normal     
x <- rlnorm(5000,meanlog = meanL,sdlog = sdL) # values    
oldpar <- parset(plots=c(1,2)) # simplifies plots par() definition    
hist(x[x < 8.0],breaks=seq(0,8,0.25),col=0,main="")     
meanx <- mean(log(x)); sdx <- sd(log(x))    
outstat <- c(exp(meanx-(sdx^2)),exp(meanx),exp(meanx+(sdx^2)/2))    
abline(v=outstat,col=c(4,1,2),lwd=3,lty=c(1,2,3))    
legend("topright",c("mode","median","bias-correct"),    
       col=c(4,1,2),lwd=3,bty="n",cex=1.2,lty=c(1,2,3))    
outh <- hist(log(x),breaks=30,col=0,main="")   # approxnormal    
hans <- addnorm(outh,log(x)) #MQMF function; try  ?addnorm    
lines(hans$x,hans$y,lwd=3,col=1) # type addnorm into the console 
par(oldpar)  # return par to old settings; this line not in book   


# R-chunk 23  Page 101  
#examine log-normal propoerties. It is a bad idea to reuse     

set.seed(12345) #'random' seeds, use getseed() for suggestions    
meanL <- 0.7;   sdL <- 0.5  #5000 random log-normal values then    
x <- rlnorm(5000,meanlog = meanL,sdlog = sdL) #try with only 500     
meanx <- mean(log(x)); sdx <- sd(log(x))    
cat("               Original  Sample \n")    
cat("Mode(x)     = ",exp(meanL - sdL^2),outstat[1],"\n")    
cat("Median(x)   = ",exp(meanL),outstat[2],"\n")    
cat("Mean(x)     = ",exp(meanL + (sdL^2)/2),outstat[3],"\n")    
cat("Mean(log(x) =  0.7     ",meanx,"\n")    
cat("sd(log(x)   =  0.5     ",sdx,"\n")    

### Fitting a Curve using Log-Normal Likelihoods    
# R-chunk 24  Page 103 
# fit a Beverton-Holt recruitment curve to tigers data Table 4.2    

data(tigers)   # use the tiger prawn data set    
lbh <- function(p,biom) return(log((p[1]*biom)/(p[2] + biom)))    
#note we are returning the log of Beverton-Holt recruitment   
pars <- c("a"=25,"b"=4.5,"sigma"=0.4)   # includes a sigma    
best <- nlm(negNLL,pars,funk=lbh,observed=log(tigers$Recruit),    
            biom=tigers$Spawn,typsize=magnitude(pars))    
outfit(best,backtran=FALSE,title="Beverton-Holt Recruitment")    
predR <- exp(lbh(best$estimate,tigers$Spawn))     
#note exp(lbh(...)) is the median because no bias adjustment    
result <- cbind(tigers,predR,tigers$Recruit/predR)    

# R-chunk 25  Page 103 
# Fig 4.12 visual examination of the fit to the tigers data    

oldp <- plot1(tigers$Spawn,predR,xlab="Spawning Biomass","Recruitment",    
      maxy=getmax(c(predR,tigers$Recruit)),lwd=2)    
points(tigers$Spawn,tigers$Recruit,pch=16,cex=1.1,col=2)  
par(oldp)  # return par to old settings; this line not in book   

# R-chunk 26  page 104, Table 4.12, code not shown in book. 
#tabulating observed, predicted and residual recruitment    

colnames(result) <- c("SpawnB","Recruit","PredR","Residual")    
kable(result,digits=c(1,1,3,4), caption='(ref:tab402)')    

### Fitting a Dynamic Model using Log-Normal Errors    
# R-chunk 27  Page 106

data(abdat)  # plot abdat fishery data using a MQMF helper  Fig 4.13    
plotspmdat(abdat) # function to quickly plot catch and cpue  


# R-chunk 28  Pages 106 - 107 
# Use log-transformed parameters for increased stability when    
# fitting the surplus production model to the abdat data-set    

param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05))     
obslog <- log(abdat$cpue) #input log-transformed observed data    
bestmod <- nlm(f=negLL,p=param,funk=simpspm,indat=as.matrix(abdat),    
               logobs=obslog)  # no typsize, or iterlim needed    
#backtransform estimates, outfit's default, as log-transformed     
outfit(bestmod,backtran = TRUE,title="abdat")        # in param    


# R-chunk 29  Pages 107 - 108 
# Fig 4.14 Examine fit of predicted to data    

predce <- simpspm(bestmod$estimate,abdat) #compare obs vs pred    
ymax <- getmax(c(predce,obslog))    
oldp <- plot1(abdat$year,obslog,type="p",maxy=ymax,ylab="Log(CPUE)",    
      xlab="Year",cex=0.9)    
lines(abdat$year,predce,lwd=2,col=2) 
par(oldp)  # return par to old settings; this line not in book    


## Likelihoods from the Binomial Distribution    
### An Example using Binomial Likelihoods    
# R-chunk 30  Page 109 
#Use Binomial distribution to test biased sex-ratio Fig 4.15    

n <- 60    # a sample of 60 animals    
p <- 0.5   # assume a sex-ration of 1:1     
m <- 1:60  # how likely is each of the 60 possibilites?    
binom <- dbinom(m,n,p)   # get individual likelihoods    
cumbin <- pbinom(m,n,p)  # get cumulative distribution    
oldp <- plot1(m,binom,type="h",xlab="Number of Males",ylab="Probability")     
abline(v=which.closest(0.025,cumbin),col=2,lwd=2) # lower 95% CI    
par(oldp)  # return par to old settings; this line not in book  


# R-chunk 31  Page 111 
# plot relative likelihood of different p values Fig 4.16    

n <- 60  # sample size; should really plot points as each independent     
m <- 20  # number of successes = finding a male    
p <- seq(0.1,0.6,0.001) #range of probability we find a male     
lik <- dbinom(m,n,p)    # R function for binomial likelihoods    
oldp <- plot1(p,lik,type="l",xlab="Prob. of 20 Males",ylab="Prob.")    
abline(v=p[which.max(lik)],col=2,lwd=2) # try "p" instead of "l" 
par(oldp)  # return par to old settings; this line not in book     

# R-chunk 32  Page 111 
# find best estimate using optimize to finely search an interval    

n <- 60; m <- 20  # trials and successes    
p <- c(0.1,0.6) #range of probability we find a male     
optimize(function(p) {dbinom(m,n,p)},interval=p,maximum=TRUE)    

### Open Bay Juvenile Fur Seal Population Size    
# R-chunk 33  Page 112 
# Juvenile furseal data-set Greaves, 1992.  Table 4.3    

furseal <- c(32,222,1020,704,1337,161.53,31,181,859,593,1125,    
             135.72,29,185,936,634,1238,153.99)    
columns <- c("tagged(m)","Sample(n)","Population(X)",    
             "95%Lower","95%Upper","StErr")    
furs <- matrix(furseal,nrow=3,ncol=6,dimnames=list(NULL,columns),    
               byrow=TRUE)    
#tabulate fur seal data Table 4.3    
kable(furs, caption='(ref:tab403)')    

# R-chunk 34  Pages 113 - 114 
# analyse two pup counts 32 from 222, and 31 from 181, rows 1-2 in    
# Table 4.3.   Now set-up storage for solutions    

optsol <- matrix(0,nrow=2,ncol=2,    
                 dimnames=list(furs[1:2,2],c("p","Likelihood")))    
X <- seq(525,1850,1) # range of potential population sizes    
p <- 151/X  #range of proportion tagged; 151 originally tagged    
m <- furs[1,1] + 1 #tags observed, with Bailey's adjustment    
n <- furs[1,2] + 1 # sample size with Bailey's adjustment    
lik1 <- dbinom(m,n,p) # individaul likelihoods    
#find best estimate with optimize to finely search an interval    
#use unlist to convert the output list into a vector    
#Note use of Bailey's adjustment (m+1), (n+1) Caughley, (1977)    
optsol[1,] <- unlist(optimize(function(p) {dbinom(m,n,p)},p,    
                              maximum=TRUE))    
m <- furs[2,1]+1;  n <- furs[2,2]+1 #repeat for sample2    
lik2 <- dbinom(m,n,p)      
totlik <- lik1 * lik2 #Joint likelihood of 2 vectors    
optsol[2,] <- unlist(optimize(function(p) {dbinom(m,n,p)},p,    
                              maximum=TRUE))    

# R-chunk 35  Page 114
# Compare outcome for 2 independent seal estimates Fig 4.17    
# Should plot points not a line as each are independent     

oldp <- plot1(X,lik1,type="l",xlab="Total Pup Numbers",    
      ylab="Probability",maxy=0.085,lwd=2)    
abline(v=X[which.max(lik1)],col=1,lwd=1)    
lines(X,lik2,lwd=2,col=2,lty=3)  # add line to plot    
abline(v=X[which.max(lik2)],col=2,lwd=1) # add optimum    
#given p = 151/X, then X = 151/p and p = optimum proportion     
legend("topright",legend=round((151/optsol[,"p"])),col=c(1,2),lwd=3,    
       bty="n",cex=1.1,lty=c(1,3)) 
par(oldp)  # return par to old settings; this line not in book     

### Using Multiple Independent Samples    

# R-chunk 36  Pages 114 - 115 
#Combined likelihood from 2 independent samples Fig 4.18    

totlik <- totlik/sum(totlik) # rescale so the total sums to one    
cumlik <- cumsum(totlik) #approx cumulative likelihood for CI        
oldp <- plot1(X,totlik,type="l",lwd=2,xlab="Total Pup Numbers",    
      ylab="Posterior Joint Probability")    
percs <- c(X[which.closest(0.025,cumlik)],X[which.max(totlik)],    
           X[which.closest(0.975,cumlik)])    
abline(v=percs,lwd=c(1,2,1),col=c(2,1,2))    
legend("topright",legend=percs,lwd=c(2,4,2),bty="n",col=c(2,1,2),    
       cex=1.2)  # now compare with averaged count    
m <- furs[3,1];  n <- furs[3,2] # likelihoods for the     
lik3 <- dbinom(m,n,p)            # average of six samples    
lik4 <- lik3/sum(lik3)  # rescale for comparison with totlik    
lines(X,lik4,lwd=2,col=3,lty=2) #add 6 sample average to plot 
par(oldp)  # return par to old settings; this line not in book     

### Analytical Approaches    
## Other Distributions    
## Likelihoods from the Multinomial Distribution    
### Using the Multinomial Distribution    
# R-chunk 37  Page 119 
#plot counts x shell-length of 2 cohorts   Figure 4.19    

cw <- 2  # 2 mm size classes, of which mids are the centers    
mids <- seq(8,54,cw) #each size class = 2 mm as in 7-9, 9-11, ...    
obs <- c(0,0,6,12,35,40,29,23,13,7,10,14,11,16,11,11,9,8,5,2,0,0,0,0)    
# data from (Helidoniotis and Haddon, 2012)    
dat <- as.matrix(cbind(mids,obs)) #xy matrix needed by inthist    
oldp <- parset()  #set up par declaration then use an MQMF function     
inthist(dat,col=2,border=3,width=1.8, #histogram of integers    
        xlabel="Shell Length mm",ylabel="Frequency",xmin=7,xmax=55)   
par(oldp)  # return par to old settings; this line not in book    

# R-chunk 38  Page 121
#cohort data with 2 guess-timated normal curves Fig 4.20    

oldp <- parset()  # set up the required par declaration    
inthist(dat,col=0,border=8,width=1.8,xlabel="Shell Length mm",    
        ylabel="Frequency",xmin=7,xmax=55,lwd=2)  # MQMF function          
#Guess normal parameters and plot those curves on histogram    
av <- c(18.0,34.5)    # the initial trial and error means and    
stdev <- c(2.75,5.75)  # their standard deviations    
prop1 <- 0.55       #  proportion of observations in cohort 1    
n <- sum(obs) #262 observations, now calculate expected counts    
cohort1 <- (n*prop1*cw)*dnorm(mids,av[1],stdev[1]) # for each    
cohort2 <- (n*(1-prop1)*cw)*dnorm(mids,av[2],stdev[2])# cohort    
#(n*prop1*cw) scales likelihoods to suit the 2mm class width    
lines(mids,cohort1,lwd=2,col=1)    
lines(mids,cohort2,lwd=2,col=4) 
par(oldp)  # return par to old settings; this line not in book     

# R-chunk 39 Page 122 
#wrapper function for calculating the multinomial log-likelihoods    
#using predfreq and mnnegLL, Use ? and examine their code    

wrapper <- function(pars,obs,sizecl,midval=TRUE) {    
  freqf <- predfreq(pars,sum(obs),sizecl=sizecl,midval=midval)    
  return(mnnegLL(obs,freqf))    
} # end of wrapper which uses MQMF::predfreq and MQMF::mnnegLL    
mids <- seq(8,54,2) # each size class = 2 mm as in 7-9, 9-11, ...    
av <- c(18.0,34.5)   # the trial and error means and    
stdev <- c(2.95,5.75)  # standard deviations    
phi1 <- 0.55      # proportion of observations in cohort 1    
pars <-c(av,stdev,phi1)  # combine parameters into a vector    
wrapper(pars,obs=obs,sizecl=mids) # calculate total -veLL    

# R-chunk 40  Page 122 
# First use the midpoints    

bestmod <- nlm(f=wrapper,p=pars,obs=obs,sizecl=mids,midval=TRUE,     
               typsize=magnitude(pars))    
outfit(bestmod,backtran=FALSE,title="Using Midpts"); cat("\n")    
#Now use the size class bounds and cumulative distribution    
#more sensitive to starting values, so use best pars from midpoints    
X <- seq((mids[1]-cw/2),(tail(mids,1)+cw/2),cw)    
bestmodb <- nlm(f=wrapper,p=bestmod$estimate,obs=obs,sizecl=X,    
                midval=FALSE,typsize=magnitude(pars))    
outfit(bestmodb,backtran=FALSE,title="Using size-class bounds")     

# R-chunk 41 Page 123 
#prepare the predicted Normal distribution curves    

pars <- bestmod$estimate # best estimate using mid-points    
cohort1 <- (n*pars[5]*cw)*dnorm(mids,pars[1],pars[3])     
cohort2 <- (n*(1-pars[5])*cw)*dnorm(mids,pars[2],pars[4])     
parsb <- bestmodb$estimate # best estimate with bounds    
nedge <- length(mids) + 1  # one extra estimate    
cump1 <- (n*pars[5])*pnorm(X,pars[1],pars[3])#no need to rescale    
cohort1b <- (cump1[2:nedge] - cump1[1:(nedge-1)])     
cump2 <- (n*(1-pars[5]))*pnorm(X,pars[2],pars[4])  # cohort 2    
cohort2b <- (cump2[2:nedge] - cump2[1:(nedge-1)])    

# R-chunk 42  Page 123 
#plot the alternate model fits to cohorts  Fig 4.21    

oldp <- parset() #set up required par declaration; then plot curves    
pick <- which(mids < 28)    
inthist(dat[pick,],col=0,border=8,width=1.8,xmin=5,xmax=28,    
        xlabel="Shell Length mm",ylabel="Frequency",lwd=3)     
lines(mids,cohort1,lwd=3,col=1,lty=2) # have used setpalette("R4")    
lines(mids,cohort1b,lwd=2,col=4)      # add the bounded results    
label <- c("midpoints","bounds")      # very minor differences   
legend("topleft",legend=label,lwd=3,col=c(1,4),bty="n",    
       cex=1.2,lty=c(2,1))    
par(oldp)  # return par to old settings; this line not in book  

# R-chunk 43  Page 124 
# setup table of results for comparison of fitting strategies    

predmid <- rowSums(cbind(cohort1,cohort2))    
predbnd <- rowSums(cbind(cohort1b,cohort2b))    
result <- as.matrix(cbind(mids,obs,predmid,predbnd,predbnd-predmid))    
colnames(result) <- c("mids","Obs","Predmid","Predbnd","Difference")    
result <- rbind(result,c(NA,colSums(result,na.rm=TRUE)[2:5]))    


# R-chunk 44  page 125, Table 4.4, code not shown in book. 
#tabulate the results of fitting cohort data  in two ways    

kable(result,digits=c(0,0,4,4,4),align=c("r","r","r","r","r"))    

## Likelihoods from the Gamma Distribution    

# R-chunk 45  Pages 126 - 127 
#Illustrate different Gamma function curves  Figure 4.22    

X <- seq(0.0,10,0.1) #now try different shapes and scale values    
dg <- dgamma(X,shape=1,scale=1)     
oldp <- plot1(X,dg,xlab = "Quantile","Probability Density")    
lines(X,dgamma(X,shape=1.5,scale=1),lwd=2,col=2,lty=2)    
lines(X,dgamma(X,shape=2,scale=1),lwd=2,col=3,lty=3)    
lines(X,dgamma(X,shape=4,scale=1),lwd=2,col=4,lty=4)    
legend("topright",legend=c("Shape 1","Shape 1.5","Shape 2",    
                           "Shape 4"),lwd=3,col=c(1,2,3,4),bty="n",cex=1.25,lty=1:4)    
mtext("Scale c = 1",side=3,outer=FALSE,line=-1.1,cex=1.0,font=7)    
par(oldp)  # return par to old settings; this line not in book  

## Likelihoods from the Beta Distribution    
# R-chunk 46  Pages 127 - 128
#Illustrate different Beta function curves. Figure 4.23    

x <- seq(0, 1, length = 1000)    
oldp <- parset()    
plot(x,dbeta(x,shape1=3,shape2=1),type="l",lwd=2,ylim=c(0,4),    
     yaxs="i",panel.first=grid(), xlab="Variable 0 - 1",     
     ylab="Beta Probability Density - Scale1 = 3")    
bval <- c(1.25,2,4,10)    
for (i in 1:length(bval))     
  lines(x,dbeta(x,shape1=3,shape2=bval[i]),lwd=2,col=(i+1),lty=c(i+1))    
legend(0.5,3.95,c(1.0,bval),col=c(1:7),lwd=2,bty="n",lty=1:5)    
par(oldp)  # return par to old settings; this line not in book  


## Bayes' Theorem    
### Introduction    
### Bayesian Methods    
### Prior Probabilities    
# R-chunk 47  Pages 132 - 133 
# can prior probabilities ever be uniniformative?  Figure 4.24    

x <- 1:1000    
y <- rep(1/1000,1000)    
cumy <- cumsum(y)    
group <- sort(rep(c(1:50),20))    
xlab <- seq(10,990,20)  
oldp <- par(no.readonly=TRUE)  # this line not in book  
par(mfrow=c(2,1),mai=c(0.45,0.3,0.05,0.05),oma=c(0.0,1.0,0.0,0.0))     
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)      
yval <- tapply(y,group,sum)    
plot(x,cumy,type="p",pch=16,cex=0.5,panel.first=grid(),    
     xlim=c(0,1000),ylim=c(0,1),ylab="",xlab="Linear Scale")    
plot(log(x),cumy,type="p",pch=16,cex=0.5,panel.first=grid(),    
     xlim=c(0,7),xlab="Logarithmic Scale",ylab="")    
mtext("Cumulative Probability",side=2,outer=TRUE,cex=0.9,font=7)
par(oldp)  # return par to old settings; this line not in book      

## End(Not run)