#' Build TMB object for geostatistical delta-GLMM
#'
#' \code{Build_TMB_Fn} builds a tagged list with everything necessary to run or interpret inputs to a geostatistical delta-GLMM
#'
#' @param TmbData, a tagged list of data inputs generated by \code{Data_Fn}
#' @param VesselConfig, a vector of form c("Vessel"=0,"VesselYear"=0) turning off vessel or vessel-year effects
#' @param Q_Config, a boolean whether to estimate catchability covariates (Q_ik in the TmbData input) or not
#' @param CovConfig, a boolean whether to estimate density covariates (X_xj and X_xtp in the TmbData input) or not
#' @param Method Spatial method used for estimation (determines bounds for logkappa)
#' @param ConvergeTol, OPTIONAL override for TMB convergence criteria
#' @param Use_REML, OPTIONAL boolean whether to use maximum marginal likelihood or restricted maximum likelihood (termed "REML")
#' @param loc_x OPTIONAL, location for each sample used to generate plausible bounds for scale parameter
#' @param Parameters OPTIONAL, a tagged list of starting parameters
#' @param Random OPTIONAL, a character vector of random effects
#' @param Map OPTIONAL, a tagged list of parameters to either mirror or turn off
#' @param DiagnosticDir OPTIONAL, a directory where diagonstic runtime information should be stored
#' @param TmbDir OPTIONAL, a directory where the CPP file for the VAST model can be found locally
#' @param RunDir OPTIONAL, a directory where the CPP file is copied, copiled, and run (must have write privileges or else the function will crash)
#' @param silent Boolean, whether TMB should run silently
#' @inheritParams Data_Fn
#' @return Tagged list containing objects for running a VAST model
#' \describe{
#' \item{Obj}{The built TMB object}
#' \item{Upper}{A vector of upper bounds for parameters, optionally for use during optimization}
#' \item{Lower}{A vector of lower bounds for parameters, optionally for use during optimization}
#' \item{Parameters}{A tagged list of parameter starting values used when building Obj, which can be extracted, modified, and then put back into \code{Build_TMB_Fn} to define different starting values}
#' \item{Map}{A taggged list of parameters to be turned off or mirrored, for similar use as Parameters}
#' \item{Random}{A character vector of random effects, for similar use as Parameters}
#' }
#' @importFrom grDevices colorRampPalette dev.off jpeg png rainbow rgb tiff
#' @importFrom graphics abline arrows axis box hist legend lines matplot mtext par plot plot.new points polygon rect title
#' @importFrom stats cor dist na.omit plogis pnorm qlogis quantile rbinom rgamma rlnorm rnorm runif sd var
#' @importFrom utils capture.output write.csv write.table
#' @export
Build_TMB_Fn <-
function( TmbData, Version, VesselConfig=c("Vessel"=0,"VesselYear"=0), Q_Config=TRUE, CovConfig=TRUE,
RhoConfig=c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsilon2"=0), Method="Mesh",
ConvergeTol=1, Use_REML=FALSE, loc_x=NULL, Parameters="generate", Random="generate", Map="generate",
DiagnosticDir=NULL, TmbDir=system.file("executables",package="SpatialDeltaGLMM"), RunDir=getwd(), silent=FALSE ){
# Compile TMB software
if( !file.exists(paste0(TmbDir,"/",Version,".cpp")) ){
stop("Make sure that ",TmbDir," contains ",Version,".cpp")
}
file.copy( from=paste0(TmbDir,"/",Version,".cpp"), to=paste0(RunDir,"/",Version,".cpp"), overwrite=FALSE)
setwd( RunDir )
compile( paste0(Version,".cpp") )
# Local functions
boundsifpresent_fn = function( par, map, name, lower, upper, bounds ){
if( name %in% names(par) ){
bounds[grep(name,names(par)),c('Lower','Upper')] = c(lower,upper)
}
return( bounds )
}
# Parameters
if( length(Parameters)==1 && Parameters=="generate" ) Parameters = Param_Fn( Version=Version, DataList=TmbData, RhoConfig=RhoConfig )
# Which parameters are turned off
if( length(Map)==1 && Map=="generate" ) Map = Make_Map( Version=Version, TmbData=TmbData, VesselConfig=VesselConfig, CovConfig=CovConfig, Q_Config=Q_Config, RhoConfig=RhoConfig, Aniso=TmbData[['Options_vec']]['Aniso'])
if( "hyperparameters_z"%in%names(Parameters) && TmbData$Options_vec['AreaAbundanceCurveTF']==0 ) Map[["hyperparameters_z"]] = factor( rep(NA,length(Parameters[["hyperparameters_z"]])) )
# Which are random
if( length(Random)==1 && Random=="generate" ){
Random = c("Epsiloninput1_st", "Omegainput1_s", "Epsiloninput2_st", "Omegainput2_s", "nu1_v", "nu2_v", "nu1_vt", "nu2_vt")
if( RhoConfig[["Beta1"]]!=0 ) Random = c(Random, "beta1_t")
if( RhoConfig[["Beta2"]]!=0 ) Random = c(Random, "beta2_t")
if( Use_REML==TRUE ){
Random = union(Random, c("beta1_t","gamma1_j","gamma1_tp","lambda1_k","beta2_t","gamma2_j","gamma2_tp","lambda2_k"))
Random = Random[which(Random %in% names(Parameters))]
}
# Avoid problems with mapping
Random = setdiff(Random, names(Map))
if( length(Random)==0) Random = NULL
}
# Build object
dyn.load( paste0(RunDir,"/",TMB::dynlib(Version)) ) # random=Random,
Obj <- MakeADFun(data=TmbData, parameters=Parameters, hessian=FALSE, map=Map, random=Random, DLL=Version, silent=silent) #
Obj$control <- list(trace=1, parscale=1, REPORT=1, reltol=1e-12, maxit=100)
# Diagnostic functions (optional)
if( !is.null(DiagnosticDir) ){
Obj$gr_orig = Obj$gr
Obj$fn_orig = Obj$fn
Obj$fn = function( vec ){
capture.output( matrix(vec,ncol=1,dimnames=list(names(Obj$par),NULL)), file=paste0(DiagnosticDir,"fn.txt") )
write.table( matrix(vec,nrow=1), row.names=FALSE, sep=",", col.names=FALSE, append=TRUE, file=paste0(DiagnosticDir,"trace.csv"))
return( Obj$fn_orig(vec) )
}
Obj$gr = function( vec ){
capture.output( matrix(vec,ncol=1,dimnames=list(names(Obj$par),NULL)), file=paste0(DiagnosticDir,"gr.txt") )
return( Obj$gr_orig(vec) )
}
write.table( matrix(Obj$par,nrow=1), row.names=FALSE, sep=",", col.names=FALSE, file=paste0(DiagnosticDir,"trace.csv"))
}
# Declare upper and lower bounds for parameter search
Bounds = matrix( NA, ncol=2, nrow=length(Obj$par), dimnames=list(names(Obj$par),c("Lower","Upper")) )
Bounds[,'Lower'] = rep(-50, length(Obj$par))
Bounds[,'Upper'] = rep( 50, length(Obj$par))
Bounds[grep("logsigmaV",names(Obj$par)),'Lower'] = log(0.01)
Bounds[grep("logsigmaVT",names(Obj$par)),'Lower'] = log(0.01)
Bounds[grep("logtau",names(Obj$par)),'Upper'] = 10 # Version < v2i
Bounds[grep("logeta",names(Obj$par)),'Upper'] = log(1/(1e-2*sqrt(4*pi))) # Version >= v2i: Lower bound on margSD = 1e-4
Bounds[grep("SigmaM",names(Obj$par)),'Upper'] = 10 # ZINB can crash if it gets > 20
if( !is.null(loc_x) && !is.na(TmbData$Options_vec['Method']) && TmbData$Options_vec['Method']==0 && Method!="Spherical_mesh" ){
Dist = dist(loc_x)
Bounds[grep("logkappa",names(Obj$par)),'Lower'] = log( sqrt(8)/max(Dist) ) # Range = nu*sqrt(8)/kappa
Bounds[grep("logkappa",names(Obj$par)),'Upper'] = log( sqrt(8)/min(Dist) ) # Range = nu*sqrt(8)/kappa
}
if( !is.na(TmbData$Options_vec['Method']) && TmbData$Options_vec['Method']==1 && Method!="Spherical_mesh" ){
Bounds[grep("logkappa",names(Obj$par)),'Upper'] = log(0.9999) # Has to be negative because exp(logkappa)=Rho < 1
}
Bounds = boundsifpresent_fn( par=Obj$par, name="gamma1", lower=-20, upper=20, bounds=Bounds)
Bounds = boundsifpresent_fn( par=Obj$par, name="gamma2", lower=-20, upper=20, bounds=Bounds)
Bounds = boundsifpresent_fn( par=Obj$par, name="lambda1", lower=-20, upper=20, bounds=Bounds)
Bounds = boundsifpresent_fn( par=Obj$par, name="lambda2", lower=-20, upper=20, bounds=Bounds)
Bounds = boundsifpresent_fn( par=Obj$par, name="Beta_rho1", lower=-0.99, upper=0.99, bounds=Bounds)
Bounds = boundsifpresent_fn( par=Obj$par, name="Beta_rho2", lower=-0.99, upper=0.99, bounds=Bounds)
Bounds = boundsifpresent_fn( par=Obj$par, name="Epsilon_rho1", lower=-0.99, upper=0.99, bounds=Bounds)
Bounds = boundsifpresent_fn( par=Obj$par, name="Epsilon_rho2", lower=-0.99, upper=0.99, bounds=Bounds)
# Change convergence tolerance
Obj$env$inner.control$step.tol <- c(1e-8,1e-12,1e-15)[ConvergeTol] # Default : 1e-8 # Change in parameters limit inner optimization
Obj$env$inner.control$tol10 <- c(1e-6,1e-8,1e-12)[ConvergeTol] # Default : 1e-3 # Change in pen.like limit inner optimization
Obj$env$inner.control$grad.tol <- c(1e-8,1e-12,1e-15)[ConvergeTol] # # Default : 1e-8 # Maximum gradient limit inner optimization
# Print number of parameters
message("Number of fixed and random effects:")
print( table(names(Obj$env$last.par)) )
# Return stuff
Return = list("Obj"=Obj, "Upper"=Bounds[,'Upper'], "Lower"=Bounds[,'Lower'], "Parameters"=Parameters, "Map"=Map, "Random"=Random)
return( Return )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.