q_model branch prior to testing. Need to merge in changes in devel

R/asap3_to_wham.R

@@ -151,7 +151,7 @@ asap3_to_wham <- function(asap3, recruit_model = 2, model_name)
   data$bias_correct_pe = 0 #bias correct log-normal process errors?
   data$bias_correct_oe = 0 #bias correct log-normal observation errors?
   data$Fbar_ages = 1:data$n_ages
-  data$simulate_state = rep(1,4) #simulate state variables (NAA, M, sel, Ecov)
+  data$simulate_state = rep(1,5) #simulate state variables (NAA, M, sel, Ecov, q)
   data$simulate_data = rep(1,3) #simulate data types (catch, indices, Ecov)
   data$simulate_period = rep(1,2) #simulate above items for (model years, projection years)
   data$percentSPR = 40 #percentage of unfished SSB/R to use for SPR-based reference points

R/prepare_projection.R

@@ -87,16 +87,17 @@ prepare_projection = function(model, proj.opts)
 
   # check ecov options are valid
   if(any(input1$data$Ecov_model > 0)){
+    n_effects = dim(input1$par$Ecov_beta)[1]
     # if Ecovs already extend beyond population model, Ecov proj options are unnecessary
-    n.beyond <- rep(NA, data$n_Ecov)
-    end.beyond <- rep(NA, data$n_Ecov)
-    for(i in 1:data$n_Ecov){
-      n.beyond[i] <- data$n_years_Ecov-1-data$ind_Ecov_out_end[i]
+    n.beyond <- end.beyond <- integer()
+    for(i in 1:data$n_Ecov) {
+      n.beyond[i] = data$n_years_Ecov-1-max(data$ind_Ecov_out_end[i,])
       end.beyond[i] <- min(n.beyond[i], proj.opts$n.yrs)     
       if(end.beyond[i] == proj.opts$n.yrs) print(paste0("ecov ",i," already fit through projection years. Using fit ecov ",i," for projections..."))
     }
     Ecov.proj <- matrix(rep(NA, max(proj.opts$n.yrs-end.beyond)*data$n_Ecov), ncol=data$n_Ecov)
-    Ecov.map <- matrix(rep(NA, max(proj.opts$n.yrs-end.beyond)*data$n_Ecov), ncol=data$n_Ecov)
+    #this isn't used anywhere
+    #Ecov.map <- matrix(rep(NA, max(proj.opts$n.yrs-end.beyond)*data$n_Ecov), ncol=data$n_Ecov)
     if(all(end.beyond < proj.opts$n.yrs)){ # if Ecov proj options ARE necessary, check they are valid
       Ecov.opt.ct <- sum(proj.opts$cont.ecov, proj.opts$use.last.ecov, !is.null(proj.opts$avg.ecov.yrs), !is.null(proj.opts$proj.ecov))
       if(Ecov.opt.ct == 0) proj.opts$cont.ecov = TRUE; Ecov.opt.ct = 1;
@@ -210,14 +211,24 @@ prepare_projection = function(model, proj.opts)
     map$log_NAA = factor(tmp)
   }
 
+    # pad q_re
+    par$q_re <- rbind(par$q_re[1:data$n_years_model,,drop=F], matrix(0, data$n_years_proj, data$n_indices))
+    map$q_re <- rbind(matrix(as.integer(map$q_re), data$n_years_model, data$n_indices), matrix(NA, data$n_years_proj, data$n_indices))
+    if(any(data$use_q_re == 1)){
+      ind = which(data$use_q_re)
+      map$q_re[,ind] <- 1:(length(ind) * (data$n_years_model + data$n_years_proj)) #turn on appropriate columns of q_re
+    }
+    map$q_re <- factor(map$q_re)
+
+
   # pad Ecov_re for projections
   if(any(data$Ecov_model > 0)){
     if(any(end.beyond < proj.opts$n.yrs)){ # need to pad Ecov_re
       for(j in 1:data$n_Ecov){
         # for(i in 1:(proj.opts$n.yrs-end.beyond[j])){
         for(i in 1:max(proj.opts$n.yrs-end.beyond)){
           if(!is.null(proj.opts$use.last.ecov)) if(proj.opts$use.last.ecov){ # use last Ecov (pad Ecov_re but map to NA)
-            Ecov.proj[i,j] <- model$rep$Ecov_re[data$ind_Ecov_out_end[j]+1+end.beyond[j],j]
+            Ecov.proj[i,j] <- model$rep$Ecov_re[max(data$ind_Ecov_out_end[j,])+1+end.beyond[j],j]
           }
           if(!is.null(proj.opts$avg.ecov.yrs)){ # use average Ecov (pad Ecov_re but map to NA)
             ecov.yrs <- data$year1_Ecov+0:(data$n_years_Ecov-1+data$n_years_proj_Ecov)

R/prepare_wham_input.R

@@ -38,12 +38,18 @@
 #'     \item{$year}{Years corresponding to observations (vector of same length as \code{$mean} and \code{$logsigma})}
 #'     \item{$use_obs}{T/F (or 0/1) vector/matrix of the same dimension as \code{$mean} and \code{$logsigma}.
 #'     Use the observation? Can be used to ignore subsets of the ecov without changing data files.}
-#'     \item{$lag}{Offset between the ecov observations and their affect on the stock.
-#'     I.e. if SST in year \emph{t} affects recruitment in year \emph{t + 1}, set \code{lag = 1}.}
+#'     \item{$lag}{integer vector of offsets between the ecov observation/process and their affect on the stock.
+#'     I.e. if SST in year \emph{t} affects recruitment in year \emph{t + 1}, set \code{lag = 1}. May also be a list (length=n_Ecov) of vectors (lenght = 2+n_indices) if multiple effects of one or more Ecovs are modeled.}
 #'     \item{$process_model}{Process model for the ecov time-series. \code{"rw"} = random walk, \code{"ar1"} = 1st order autoregressive, \code{NA} = do not fit}
-#'     \item{$where}{Where does the ecov affect the population? \code{"recuit"} = recruitment,
-#'     \code{"M"} = natural mortality, \code{"none"} = fit ecov but without an effect on the population.}
-#'     \item{$how}{How does the ecov affect the \code{$where} process? These options are
+#'     \item{$where}{Character vector for where each ecov affects the population. \code{"recuit"} = recruitment,
+#'     \code{"M"} = natural mortality, \code{"q"} = catchability for indices, \code{"none"} = fit ecov but without an effect on the 
+#'     population. May also be a list (element for each ecov) of multiple character vectors so each ecov can multiple effects.}
+#'     \item{$indices}{indices that each ecov affects. Must be a list of length n.ecov, where each element is a vector of indices (1:n_indices). Must be provided when any of \code{where} = "q"}
+#'     \item{$link_model}{vector or (orthogonal polynomial order) models for effects or each ecov on the \code{$where} process. Options: 'linear' (default) or 'poly-x'
+#'     where x = 2, ... (e.g. 'poly-2' specifies a quadratic model, \eqn{b_0 + b_1*ecov + b_2*ecov^2 + ...}). A list (length = n_Ecov) of character vectors for modeling
+#'      effects on recruitment, M, and catchabilities for index 1,..., index n_indices (length of the vector is 2 + n_indices).}
+#'     \item{$ages}{Ages that each ecov affects. Must be a list of length n.ecov, where each element is a vector of ages. Default = list(1:n.ages) to affect all ages.}
+#'     \item{$how}{How does the ecov affect the \code{$where} process? Currently this only isThese options are
 #'     specific to the \code{$where} process.
 #' Options for recruitment are described in \href{https://www.sciencedirect.com/science/article/pii/S1385110197000221}{Iles & Beverton (1998)}:
 #'   \describe{
@@ -58,10 +64,12 @@
 #'   \describe{
 #'     \item{= 0}{none (but fit ecov time-series model in order to compare AIC)}
 #'     \item{= 1}{effect on mean M (shared across ages)}
+#'   }
+#' Options for catchability:
+#'   \describe{
+#'     \item{= 0}{none (but fit ecov time-series model in order to compare AIC)}
+#'     \item{= 1}{effect on one or more catchabilities (see \code{$indices)})}
 #'   }}
-#'     \item{$link_model}{Model describing ecov effect on the \code{$where} process. Options: 'linear' (default) or 'poly-x'
-#'     where x = 2, ... (e.g. 'poly-2' specifies a quadratic model, \eqn{b0 + b1*ecov + b2*ecov^2 + ...}).}
-#'     \item{$ages}{Ages that each ecov affects. Must be a list of length n.ecov, where each element is a vector of ages. Default = list(1:n.ages) to affect all ages.}
 #'   }
 #'
 #' \code{selectivity} specifies options for selectivity, to overwrite existing options specified in the ASAP data file.
@@ -166,23 +174,25 @@
 #'		 }
 #'   }
 #'
-#' \code{catchability} specifies options for catchability. If \code{NULL} and \code{asap3} is not NULL, a single catchability parameter used with initial values specified in ASAP file. If both are NULL, initial catchabilities for all indices = 0.3.
+#' \code{catchability} specifies options for catchability. If \code{NULL} and \code{asap3} is not NULL, a single catchability parameter for each index is used with initial values specified in ASAP file. If both are NULL, initial catchabilities for all indices = 0.3.
 #' Otherwise, it is a list with the following optional entries:
 #'   \describe{
 #'     \item{$re}{Time-varying (random effects) for each index. Vector with length = number of indices.
-#'                  If \code{NULL}, catchability parameters are constant over time.
 #'                  Each entry of \code{catchability$re} must be one of the following options:
 #'                  \describe{
 #'                    \item{"none"}{(default) are constant}
 #'                    \item{"iid"}{vary by year and age/par, but uncorrelated}
 #'                    \item{"ar1"}{correlated by year (AR1)}
 #'                  }
 #'                 }
-#'     \item{$initial_q}{Initial catchabilities for each index. vector length = number of indices. Will override values provided in \code{basic_info$q}. If \code{basic_info$q} and \code{asap3} are not provided default q values are 0.3.}
+#'     \item{$initial_q}{Initial catchabilities for each index. vector length = number of indices. Will override values provided in \code{basic_info$q}. 
+#'        If \code{basic_info$q} and \code{asap3} are not provided, default q values are 0.3.}
 #'     \item{$q_lower}{Lower bound for catchabilities for each index. vector length = number of indices. For indices with NULL components default lower values are 0.}
 #'     \item{$q_upper}{Upper bound for catchabilities for each index. vector length = number of indices. For indices with NULL components default lower values are 1000.}
 #'     \item{$fix_q}{Which catcabilities not estimated. vector of integers indicating which indices have catchability fixed. E.g. c(2,3) when the second and third catchabilties are fixed.}
-#'     \item{$prior}{vector of NULL and sigmas to use for gaussian prior on logit transform of catchability parameter. Length = number of indices. Indices with non-NULL values will have mean logit q as a random effect with mean determined by logit transform of \code{catchability$initial_q} and sigma as standard error.}
+#'     \item{$prior_sd}{vector of NA and standard deviations to use for gaussian prior on logit transform of catchability parameter. Length = number of indices. 
+#'       Indices with non-NA values will have mean logit q as a random effect with mean determined by logit transform of \code{catchability$initial_q} and  
+#'       sigma as standard error.}
 #'   }
 #'
 #' \code{age_comp} specifies the age composition models for fleet(s) and indices.
@@ -240,7 +250,7 @@
 #'     \item{$percentFMSY}{(0-100) percentage of Fmsy to use in projections.}
 #'     \item{$XSPR_R_avg_yrs}{which years to average recruitments for calculating SPR-based SSB reference points.}
 #'     \item{$XSPR_R_opt}{1(3): use annual R estimates(predictions) for annual SSB_XSPR, 2(4): use average R estimates(predictions).}
-#'     \item{$simulate_process_error}{T/F vector (length = 4). When simulating from the model, whether to simulate any process errors for NAA, M, selectivity, and Ecov. Only used if applicable.}
+#'     \item{$simulate_process_error}{T/F vector (length = 5). When simulating from the model, whether to simulate any process errors for NAA, M, selectivity, Ecov, and q. Only used if applicable.}
 #'     \item{$simulate_observation_error}{T/F vector (length = 3). When simulating from the model, whether to simulate  catch, index, and ecov observations.}
 #'     \item{$simulate_period}{T/F vector (length = 2). When simulating from the model, whether to simulate base period (model years) and projection period.}
 #'     \item{$bias_correct_process}{T/F. Perform bias correction of log-normal random effects for NAA.}
@@ -330,7 +340,7 @@ prepare_wham_input <- function(asap3 = NULL, model_name="WHAM for unnamed stock"
 	if(is.null(catchability)){
 		q_opts = NULL
 		q_names = c("q")
-		if(any(names(basic_info) %in% q_names)) q_opts$initial_q = basic_info[q_names]
+		if(any(names(basic_info) %in% q_names)) q_opts$initial_q = basic_info[[q_names]]
 	} else {
 		q_opts = catchability
 	}
@@ -476,7 +486,7 @@ initial_input_fn = function(input, basic_info){
 
   input$data$bias_correct_pe = 1 #bias correct log-normal process errors?
   input$data$bias_correct_oe = 1 #bias correct log-normal observation errors?
-  input$data$simulate_state = rep(1,4) #simulate state variables (NAA, M, sel, Ecov)
+  input$data$simulate_state = rep(1,5) #simulate state variables (NAA, M, sel, Ecov)
   input$data$simulate_data = rep(1,3) #simulate data types (catch, indices, Ecov)
   input$data$simulate_period = rep(1,2) #simulate above items for (model years, projection years)
   input$data$percentSPR = 40 #percentage of unfished SSB/R to use for SPR-based reference points