/multi_SAM

Connecting single-stock assessment models through correlated survival

Primary LanguageR

multiStockassessment

The multiStockassessment package extends the stockassessment to allow linked stocks (See e.g. Albertsen, C. M., Nielsen, A., and Thygesen, U. H. (2018) Connecting single-stock assessment models through correlated survival. ICES Journal of Marine Science, 75(1), 235-244. doi: 10.1093/icesjms/fsx114). The package can be installed with devtools:

devtools::install_github("calbertsen/multi_SAM", subdir = "multiStockassessment")

Short example

To use the multiStockassessment package, the individual stocks must first be fitted with the stockassessment package and combined to a samset using the c(...) function.

Preparing data using stockassessment

The code below loads the North Sea cod data from the stockassessment package, fits it, simulates two new data sets and fit those individually.

library(stockassessment)

data(nscodData)
data(nscodConf)
data(nscodParameters)

set.seed(9876)
fit <- sam.fit(nscodData, nscodConf, nscodParameters,sim.condRE = FALSE)

datSim <- simulate(fit,nsim=2)
fitSim <- do.call("c",lapply(datSim,function(x)sam.fit(x,nscodConf, nscodParameters)))

Creating correlation structure

library(multiStockassessment)

After the simulated data have been fitted and combined by the c function, the multiStockassessment package can be used to fit a model with correlated survival.

The suggestCorStructure function can be used to set up the correlation matrix. The nAgeClose argument can be used to construct the banded structures used in the paper. The result is a logical symmetric matrix indicating whether a correlation should be fixed at zero. The dimension of the matrix is the number of stocks times the number of ages in the data.

For instance, suggestCorStructure(fitSim,nAgeClose=1) creates a matrix where ages less than 1 appart are correlated between the stocks:

cs <- suggestCorStructure(fitSim,nAgeClose=1)
cs
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
##  [1,]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
##  [2,]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
##  [3,]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE
##  [4,]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
##  [5,]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
##  [6,]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [7,] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [8,]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [9,]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [10,]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [11,]  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [12,]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
##       [,12]
##  [1,]  TRUE
##  [2,]  TRUE
##  [3,]  TRUE
##  [4,]  TRUE
##  [5,]  TRUE
##  [6,] FALSE
##  [7,]  TRUE
##  [8,]  TRUE
##  [9,]  TRUE
## [10,]  TRUE
## [11,]  TRUE
## [12,]  TRUE

The suggestCorStructure function has several options to configure the correlation structure, and the result can be modified by hand for complete freedom.

Fitting the model

The model is fitted by the multisam.fit function. The function requires a samset and a correlaiton structure. By default, the correlation structure given models the partial correlations between cohorts. This can be changed to regular correlations by setting usePartialCors = FALSE.

obj <- multisam.fit(fitSim,cs)
obj
## Multi-SAM model with 2 stocks: log likelihood is -253.2214. Convergence OK.

Investigating the result

Most plotting and table functions from the stockassessment package have a corresponding version in multiStockassessment.

The fitted correlation structure can be plotted by

corplot(obj)

Other plotting functions implemented from the stockassessment package are

fbarplot(obj)

ssbplot(obj)

ssbplot(obj)

tsbplot(obj)

recplot(obj)

catchplot(obj)

srplot(obj)

fitplot(obj,1)

Note that the fitplot function above requires the stock to be specified.

fitplot(obj,2)

Comparing models

Models can be compared with the modeltable function. For instance, the model above can be compared to individual assessments by

cs2 <- suggestCorStructure(fitSim,nAgeClose=0)
obj2 <- multisam.fit(fitSim,cs2)
modeltable(obj,obj2)
##       log(L) #par      AIC Pval( M1 -> M2 )
## M1 -253.2214   74 654.4428               NA
## M2 -255.0715   68 646.1429        0.7171872