diagnostics are important for replicability, by ensuring that a global solution has actually been found and that assumptions arnt violated, so when the assessment is repeated you get a similar result.
Residuals
Patterns in residuals from the fits of the CPUE to stock abundance may indicate a violation of models assumptions. Which may result in biased estimates of parameters, reference points and stock trends. While variance estimates obtained from bootstrapping assume that residuals are Independently and Identically Distributed (i.i.d.).
The residuals are in the diags slot.
print(head(bdHat@diags),digits=3)
name year obs hat residual residualLag qqx qqy qqHat
1 1 1 1061 984 0.075 -0.144 0.184 0.075 0.0279
2 1 2 846 977 -0.144 -0.181 -0.579 -0.144 -0.1514
3 1 3 800 960 -0.181 0.113 -0.708 -0.181 -0.1816
4 1 4 1051 938 0.113 -0.061 0.518 0.113 0.1066
5 1 5 861 915 -0.061 0.102 -0.237 -0.061 -0.0710
6 1 6 988 892 0.102 0.166 0.346 0.102 0.0661
Checking the distribution of residuals can be done by plotting the obsevered quantiles against the predicted quantiles from the assumed distribution using Q-Q plots. These compare a sample of data (the residuals) on the vertical axis to a statistical population (e.g. from a normal distribution) on the horizontal axis. Any nonlinear patterns may imply that the data are not normally distributed i.e. \(X ~ N(0,1)\), for example a systematic departure from a straight line may indicate skewness or over or under dispersion.
rsdl=bdHat@diags
ggplot(rsdl) +
geom_point( aes(qqx,qqy)) +
stat_smooth(aes(qqx,qqHat),method="lm",se=T,fill="blue", alpha=0.1) +
theme_bw()+theme(legend.position="bottom")
It is assumed that an index is proportional to the stock so when plotting the observed against the fitted values the points should fall around the \(y=x\) line, if they do not then the index may not be a good proxy for the stock trend.
library(diags)
ggplot(with(rsdl, data.frame(obs=stdz(obs),hat=stdz(hat)))) +
geom_abline(aes(slope=1,intercept=0)) +
geom_point( aes(obs,hat)) +
stat_smooth(aes(obs,hat),method="lm", se=F) +
theme_bw()+theme(legend.position="bottom") +
xlab("Fitted") + ylab("Observed")
To look for systematic patterns the residuals can be plotted by year, a lowess smoother helps to identify if the proxy doesnt agree with the estimated stock trend based on the catch
dat=transform(subset(rsdl,!is.na(residual),
residual=stdz(residual,na.rm=T)))
ggplot(aes(year,residual),data=dat) +
geom_hline(aes(yintercept=0)) +
geom_point() +
stat_smooth(method="loess",se=F) +
theme_bw()+theme(legend.position="bottom")
It is also assumed that variance of the index does not vary with the mean, this can be checked by plotting the residuals against the fitted values.
ggplot(aes(hat, residual),
data=subset(rsdl,!is.na(hat) & !is.na(residual))) +
geom_hline(aes(yintercept=0)) +
geom_point() +
stat_smooth(method="loess",se=F) +
theme_bw()+theme(legend.position="bottom")
It is assumed that the residuals are not autocorrelated, which can be checked by plotting the residuals against each other with a lag of 1. Significant autocorrelations could be due to an increase in catchability with time, which may result in a more optimistic estimate of current stock status as a decline in the stock may be masked by an increase in catchability.
[1] 1.92
ggplot(rsdl) +
geom_point( aes(residual,residualLag)) +
stat_smooth(aes(residual,residualLag),method="lm",se=F) +
geom_hline(aes(yintercept=0)) +
xlab(expression(Residual[t])) +
ylab(expression(Residual[t+1])) +
theme_bw()+theme(legend.position="bottom")
Warning: Removed 1 rows containing non-finite values (stat_smooth).
Warning: Removed 1 rows containing missing values (geom_point).
Profiles
Likelihood profiles are useful to check that you are actually at a global solution and not stuck on a small hill with your back to the mountain. They are also useful for evaluating the infomation content of the data and whether different data sets are telling you different things and you need to ask more questions to determine the truth.
The control slot can be used to produce a profile, i.e. fix a parameter or parameters for a range of values and then find the maximum likelihood by estimating the other parameters.
1D
bdHat=fit(bdHat,cpue)
setControl(bdHat)=params(bdHat)
res=profile(bdHat,which='r',fixed=c('b0','p'),
cpue,range=seq(0.95,1.03,.002))
ggplot(subset(res,ll<0))+
geom_line(aes(r,ll)) +
theme_bw()
res=profile(bdHat,which=c('r','k'),fixed=c('b0','p'),
cpue,range=seq(0.97,1.03,.02))
ggplot(res, aes(r, k, z=ll))+
stat_contour(aes(colour = ..level..), size = 1)+
theme_bw()
likelihood components
bd=sim()
set.seed(7890)
Us =FLQuants("Unbiased" =
rlnorm(1,log((stock(bd)[,-dims(bd)$year]+
stock(bd)[,-1])/2),0.2),
"Increase in q"=
rlnorm(1,log((stock(bd)[,-dims(bd)$year]+
stock(bd)[,-1])/2),0.2))
setParams( bd)=Us
setControl(bd)=params(bd)
bd@control[3:4,"phase"]=-1
bd=fit(bd,index=Us)
bd@control[,c("min")]=bd@params*0.1
bd@control[,c("val")]=bd@params
bd@control[,c("max")]=bd@params*10
prfl=profile(bd,which='r',index=Us,
range=seq(0.975,1.05,.001))
ggplot(prfl)+
geom_path(aes(r,ll,group=index,col=index))+
facet_wrap(~index,scale="free",ncol=1) +
theme(legend.position="bottom")+
theme_bw()
Profile Slot
\newpage
