vector_combine: combine vectors using Conn 2010 hierarchical method

In nikolaifish/bamExtras: Functions and Data in Support of Stock Assessments with the Beaufort Assessment Model (BAM)

combine vectors using Conn 2010 hierarchical method

Description

This is a version of a function that was written by Paul Conn, Kyle Shertzer, and Erik Williams, using jags to combine indices of abundance. Nikolai Klibansky converted that script into this function.

Usage

vector_combine(
  x,
  cv,
  key_nm = NULL,
  jags_args = list(parameters.to.save = c("chi", "mu", "sigma"), n.chains = 3, n.thin =
    1, n.burnin = 10000, n.iter = 80000),
  cv.na = 0.05,
  run_mcmcplot = FALSE,
  mcmcplot_args = list(dir = "mcmcplot")
)

Arguments

x	observed vectors as a data frame
cv	observed cvs associated with x, as a data frame
key_nm	If a common column, like "year", is provided in x and cv as a key to relate values between data frames provide the name as key_nm. Should be coercible to a numeric vector. If NULL, then row names will be used.
jags_args	list of arguments to be passed to R2jags::jags() within the function
cv.na	value to to replace missing cv values with
run_mcmcplot	Should mcmcplots::mcmcplot be run to plot jags results?
mcmcplot_args	if run_mcmcplot = TRUE, additional arguements can be passed to mcmcplot

Details

This function was designed to be used to combine multiple indices of abundance into a single index, but it could also be used to combine other types of numeric vectors.

Author(s)

Nikolai Klibansky, Kyle Shertzer, Erik Williams

References

Conn, P. B. 2010. Hierarchical analysis of multiple noisy abundance indices. Canadian Journal of Fisheries and Aquatic Sciences 67:108-120.

Examples

## Not run: 
# Example usage
library(ggplot2)
library(tidyr)

# Combining commercial indices of abundance from the recent blueline tilefish assessment for the US Southeast Atlantic
ts <- rdat_BluelineTilefish$t.series

x <- ts[,grepl("^U.c.*.ob$",names(ts))]
cv <- ts[,grepl("^cv.U.c",names(ts))]

names(x) <- gsub(".ob$","",names(x))
names(cv) <- gsub(".ob$","",names(cv))

# Only include years where at least one index had real values
yrs <- rownames(x)[apply(x,1,function(y){any(!is.na(y))})]
x <- x[yrs,]
cv <- cv[yrs,]

# run the function
vc_out <- vector_combine(x=x,cv=cv)

# plot results with mcmcplot (use dir argument if you want to save the plots somewhere)
mcmcplots::mcmcplot(vc_out$mcmcout)

# plot original indices and resulting index (with ggplot, to be hip)
dfout <- vc_out$dfout
x2 <- cbind(x,"combined"=dfout$mu)
sd <- x*cv
sd2 <- cbind(cv,"combined"=dfout$sd)
names(sd2) <- gsub("^cv.","",names(sd2))
x2lo <- x2-2*sd2
names(x2lo) <- paste0(names(x2lo),".lo")
x2up <- x2+2*sd2
names(x2up) <- paste0(names(x2up),".up")

x2 <- cbind("year"=rownames(x2),x2)
sd2 <- cbind("year"=rownames(sd2),sd2)

x3 <- tidyr::pivot_longer(x2,cols=!year,names_to="fleet",values_to="index")
sd3 <- tidyr::pivot_longer(sd2,cols=!year,names_to="fleet",values_to="sd")

x4 <- na.omit(merge(x3,sd3,by=c("year","fleet")))

p <- ggplot(x4, aes(x=year, y=index, group=fleet, color=fleet,linetype=fleet,size=fleet)) +
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=index-2*sd, ymax=index+2*sd), width=1, position=position_dodge(0.2)) +
  labs(title="indices of abundance", x="year", y = "index of abundance") +
  theme_classic() +
  scale_color_manual(values=c("black","red","blue")) +
  scale_linetype_manual(values=c("solid","dashed","dashed")) +
  scale_size_manual(values=c(1,0.5,0.5))

p

## End(Not run)

nikolaifish/bamExtras documentation built on July 21, 2023, 8:26 a.m.