forked from calbertsen/multi_SAM
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.Rmd
132 lines (84 loc) · 3.47 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
# multiStockassessment
The `multiStockassessment` package extends the [`stockassessment`](https://github.com/fishfollower/SAM) 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](https://dx.doi.org/10.1093/icesjms/fsx114)). The package can be installed with `devtools`:
```r
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.
```{r,results='hide',warning=FALSE}
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
```{r}
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:
```{r}
cs <- suggestCorStructure(fitSim,nAgeClose=1)
cs
```
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`.
```{r,results="hide",warning=FALSE}
obj <- multisam.fit(fitSim,cs)
```
```{r}
obj
```
### 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
```{r}
corplot(obj)
```
Other plotting functions implemented from the `stockassessment` package are
```{r}
fbarplot(obj)
```
```{r}
ssbplot(obj)
```
```{r}
ssbplot(obj)
```
```{r}
tsbplot(obj)
```
```{r}
recplot(obj)
```
```{r}
catchplot(obj)
```
```{r}
srplot(obj)
```
```{r}
fitplot(obj,1)
```
Note that the `fitplot` function above requires the stock to be specified.
```{r}
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
```{r,results="hide",warning=FALSE}
cs2 <- suggestCorStructure(fitSim,nAgeClose=0)
obj2 <- multisam.fit(fitSim,cs2)
```
```{r}
modeltable(obj,obj2)
```