| chapter4 {MQMF} | R Documentation |
chapter4 The 47 R-code chunks from Model Parameter Estimation
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
verbose |
Should instructions to written to the console, default=TRUE |
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)
[Package MQMF version 0.1.5 Index]