#Controlling R code for parent-offspring CKMR estimation #author V Trenkel #last modifed February 2020 library(TMB) #use TMB for creating derivatives of model variables #compile model compile("CKMRPOP_RanEff.cpp") dyn.load(dynlib(ModelName)) #load data load(paste("PopDataCKMR_offshore.Rdata",sep="")) #create parameter list and initialise parameters<-list( n0_parvec=20, #log(N0) numbers in year0 same number for each sex [female, male] logSdGamma=-3, #log SD for random growth parameter (random walk) gammainit=0, U=rep(0,(Data$last_y))) # random effects vector for population growth rates single value for both sexes #estimate all parameters obj <- MakeADFun(Data, parameters, DLL=ModelName,hessian=TRUE,random=c("U")) opt <- nlminb(obj$par,obj$fn,obj$gr) est <- sdreport(obj) # adrep <- summary(est,"report") #gets estimates for parameter values and reported variables ADREPORT(xx) vars=rownames(adrep) adrep=data.frame(adrep,Variable=vars) #adrep contains abundance estimates, first females, then males and negative loglikelihood id=grep("tot_lglk",rownames(adrep)) logLL=-adrep$Estimate[id] #loglikelihood id=grep("N_y",rownames(adrep)) Nest=adrep$Estimate[id] #estimated abundance SD=adrep$"Std..Error"[id] npar=length(opt$par) AIC=2*npar-2*logLL #AIC #create matrix of abundance estimates, standard deviation, 95% confidence intervals, CV Years=c(min(birthyears)-1,birthyears) Nest=data.frame(Years=Years,Nest=Nest,SD=SD,NCI2.5=Nest-1.96*SD,NCI97.5=Nest+1.96*SD,CV=SD/Nest) #plot results birthyears=Nest$Years #select estimates after 2010 id=Nest$Years>=2011 ylimy=c(0,max(Nest[id,c("NCI97.5")]))/1000 plot(birthyears[id],Nest[id,"Nest"]/1000,ylim=ylimy,type="l",col="grey",xlab="Year",ylab="Abundance of spawning individuals (thousands)",cex.lab=1, las=1,main=subdat) polygon(c(birthyears[id],rev(birthyears[id])),c(Nest[id,c("NCI2.5")],rev(Nest[id,c("NCI97.5")]))/1000, xlab="",ylab="",lty=1, col=rgb(25,25,25, max = 255, alpha = 10, names = "grey50"),ylim=ylimy,border=NA) lines(birthyears[id],Nest[id,"Nest"]/1000,col="grey",lwd=2)