andrew-edwards / EDMsimulate

An R package for simulating fish populations for empirical dynamic modeling
1 stars 1 forks source link

Should there be a check in `mve()` to make sure that the `response` argument with 0 lag is not included in the `lags` #16

Open andrew-edwards opened 11 months ago

andrew-edwards commented 11 months ago

@luke-a-rogers if you have a moment:

So the function create_subset_lags() creates all the potential lags with:

subset_lags_multi <- create_subset_lags(lags_use_multi)
subset_lags_multi
#> [[1]]
#> [[1]]$R_prime_t
#> [1] 0
#> 

etc. But we're going to have response = R_prime_t in the function call, so presumably we should not include R_prime_t values with lag of 0 when trying to estimate R_prime_t values (with lag of 0, by definition).

But we do still want to look at lags of R_prime_t, so we should use something like

lags_use_multi <- list(R_prime_t = 1:2,
                             S_t = 0:3)

with no R_prime_t = 0. I can add in a check, just checking I'm understanding it all. Think I am!

andrew-edwards commented 11 months ago

[Andy working through his answers!] Now think the above is actually wrong!

First, after talking with Carrie: We will use R_t (not R_prime_t, though will have a switch) with lags of 1, 2, and 3. S_t with lags 1 to 7. When we need the forecast for R_t we don't yet know S_t because that is calculated at the end of the year.

So, for our application, we want to exclude S_t (that's just a matter of the in order in which things happen in a year - if really wanted an S_t could define S_t as the spawners in year t-1 but that would utterly confusing and poor notation!; just need a note in the methods somewhere).

But we can include R_t (i.e. lag of 0) as we can totally use values of it (as a dimension in the state space) to look for nearest neighbours and make a projection. For focal time t*, the value of R_t* is automatically ignored as a nearest neighbour (can't be a nearest neighbour of yourself), so we are fine.

This explains why, in mve_understanding.Rmd vignette I have (currently have, may always change as I update things) observed and forecast not matching for the ssr with just R_t - of course they shouldn't, as this is just a one-dimensional state space, like left-hand part of Figure 1(c) in manuscript 1, which is not a great represenation of the dynamics. Results:

> lags_use_multi <- list(R_t = 0,
                       S_t = 0:3)
> subset_lags_multi <- create_subset_lags(lags_use_multi)
> subset_lags_multi[[1]]
$R_t
[1] 0

> mve_on_simulated$raw_forecasts[[1]] %>% as.data.frame()
   set time points dim     observed     forecast        mre      rmse R_t_0
1    0    1      0   2 1.088943e+00           NA         NA        NA     1
2    0    2      0   2 2.629785e+00           NA         NA        NA     1
3    0    3      0   2 1.459091e+00           NA         NA        NA     1
4    0    4      0   2 1.100923e+00           NA         NA        NA     1
5    0    5      0   2 1.459139e-01           NA         NA        NA     1
6    0    6      0   2 4.835006e+00           NA         NA        NA     1
7    0    7      0   2 5.302259e+00           NA         NA        NA     1
8    0    8      0   2 2.301098e+00           NA         NA        NA     1
9    0    9      0   2 2.470542e-01           NA         NA        NA     1
10   0   10      0   2 3.062190e+00           NA         NA        NA     1
11   0   11      0   2 1.152963e+00           NA         NA        NA     1
12   0   12      0   2 1.516400e+00           NA         NA        NA     1
13   0   13      0   2 2.354557e-01           NA         NA        NA     1
14   0   14      0   2 7.893669e-01           NA         NA        NA     1
15   0   15      0   2 6.122476e-01           NA         NA        NA     1
16   0   16      0   2 1.283056e+00           NA         NA        NA     1
17   0   17      0   2 3.010416e-01           NA         NA        NA     1
18   0   18      0   2 4.105536e+00           NA         NA        NA     1
19   0   19      0   2 5.162885e-01           NA         NA        NA     1
20   0   20      0   2 1.471384e+00           NA         NA        NA     1
21   0   21      0   2 1.018973e+00           NA         NA        NA     1
22   0   22      0   2 6.498114e+00           NA         NA        NA     1
23   0   23      0   2 6.903729e-01           NA         NA        NA     1
24   0   24      0   2 1.941716e+00           NA         NA        NA     1
25   0   25      0   2 2.886900e+00           NA         NA        NA     1
26   0   26      0   2 4.511924e+00           NA         NA        NA     1
27   0   27      0   2 1.309928e+00           NA         NA        NA     1
28   0   28      0   2 3.077552e-01           NA         NA        NA     1
29   0   29      0   2 2.655253e-01           NA         NA        NA     1
30   0   30      0   2 4.761460e+00           NA         NA        NA     1
31   0   31      0   2 5.050997e-01           NA         NA        NA     1
32   0   32      0   2 6.186311e-01           NA         NA        NA     1
33   0   33      0   2 2.191579e-01           NA         NA        NA     1
34   0   34      0   2 1.602444e+00           NA         NA        NA     1
35   0   35      0   2 1.722314e+00           NA         NA        NA     1
36   0   36      0   2 1.684427e+00           NA         NA        NA     1
37   0   37      0   2 6.807179e-01           NA         NA        NA     1
38   0   38      0   2 1.266255e+00           NA         NA        NA     1
39   0   39      0   2 1.669035e+00           NA         NA        NA     1
40   0   40      0   2 3.287825e+00           NA         NA        NA     1
41   0   41      0   2 3.540660e-01           NA         NA        NA     1
42   0   42      0   2 2.215093e+00           NA         NA        NA     1
43   0   43      0   2 9.327816e+00           NA         NA        NA     1
44   0   44      0   2 1.636207e+00           NA         NA        NA     1
45   0   45      0   2 3.093250e-01           NA         NA        NA     1
46   0   46      0   2 9.218332e-01           NA         NA        NA     1
47   0   47      0   2 1.255624e+01           NA         NA        NA     1
48   0   48      0   2 1.350839e+00           NA         NA        NA     1
49   0   49      0   2 2.008082e-01           NA         NA        NA     1
50   0   50      0   2 3.660299e-01           NA         NA        NA     1
51   0   51     49   2 1.012231e+00 2.146780e+00  1.1345494  1.134549     1
52   0   52     50   2 1.928335e-01 6.497960e+00  3.7198382  4.530002     1
53   0   53     51   2 7.400364e-02 4.822170e-01  2.6159632  3.706232     1
54   0   54     52   2 6.652207e-02 2.538767e+00  2.5800337  3.439494     1
55   0   55     53   2 2.199365e+00 6.684127e-02  1.6375223  3.220811     1
56   0   56     54   2 3.140110e-01 9.289672e+00  2.8605454  4.698060     1
57   0   57     55   2 2.438306e-01 9.611046e-01  2.5543638  4.357999     1
58   0   58     56   2 8.208511e-02 2.685348e+00  2.5604761  4.179146     1
59   0   59     57   2 4.759562e+00 6.747965e-01  1.8221160  4.168764     1
60   1   60     58   2 1.955115e+00 5.050997e-01  1.4949029  3.981330     1
61   1   61     59   2 6.049381e-01 2.886900e+00  1.5664537  3.857899     1
62   1   62     60   2 3.776866e-01 9.696860e-01  1.4852492  3.697609     1
63   1   63     61   2 4.171342e+01 1.320395e+00 -1.7361562 11.752786     1
64   1   64     62   2 6.872105e+00 1.247020e+00 -2.0139369 11.424615     1
65   1   65     63   2 5.456699e-01 8.042309e-01 -1.8624370 11.037428     1
66   1   66     64   2 3.576814e-01 9.930365e-01 -1.7063250 10.688124     1
67   1   67     65   2 5.051707e-02 2.292826e+00 -1.4740524 10.383256     1
68   1   68     66   2 8.461405e-03 2.015695e+00 -1.2806477 10.101796     1
69   1   69     67   2 1.369006e-04 6.870689e-01 -1.1770908  9.833629     1
70   1   70     68   2 7.435849e-05 2.231672e-03 -1.1181284  9.584636     1
71   1   71     69   2 7.879385e-03 7.435849e-05 -1.0652559  9.353647     1
72   1   72     70   2 5.787141e-02 1.369319e-04 -1.0194595  9.138600     1
73   1   73     71   2 3.647626e-03 8.744025e-01 -0.9372762  8.939571     1
74   1   74     72   2 1.563364e-04 1.964525e-02 -0.8974110  8.751349     1
75   1   75     73   2 9.082222e-03 3.748380e-04 -0.8618629  8.574536     1
76   1   76     74   2 3.212838e-02 1.637866e-02 -0.8293201  8.408025     1
77   1   77     75   2 5.279232e-03 1.325975e-02 -0.7983089  8.250852     1
78   1   78     76   2 8.900211e-04 1.659317e-02 -0.7692371  8.102177     1
79   1   79     77   2 5.969987e-03 5.648620e-03 -0.7427227  7.961258     1
80   1   80     78   2 1.096713e-01 8.549143e-03 -0.7213360  7.827468     1
   S_t_0 S_t_1 S_t_2 S_t_3
1      0     0     0     0
2      0     0     0     0
3      0     0     0     0
4      0     0     0     0
5      0     0     0     0
6      0     0     0     0
7      0     0     0     0
8      0     0     0     0
9      0     0     0     0
10     0     0     0     0
11     0     0     0     0
12     0     0     0     0
13     0     0     0     0
14     0     0     0     0
15     0     0     0     0
16     0     0     0     0
17     0     0     0     0
18     0     0     0     0
19     0     0     0     0
20     0     0     0     0
21     0     0     0     0
22     0     0     0     0
23     0     0     0     0
24     0     0     0     0
25     0     0     0     0
26     0     0     0     0
27     0     0     0     0
28     0     0     0     0
29     0     0     0     0
30     0     0     0     0
31     0     0     0     0
32     0     0     0     0
33     0     0     0     0
34     0     0     0     0
35     0     0     0     0
36     0     0     0     0
37     0     0     0     0
38     0     0     0     0
39     0     0     0     0
40     0     0     0     0
41     0     0     0     0
42     0     0     0     0
43     0     0     0     0
44     0     0     0     0
45     0     0     0     0
46     0     0     0     0
47     0     0     0     0
48     0     0     0     0
49     0     0     0     0
50     0     0     0     0
51     0     0     0     0
52     0     0     0     0
53     0     0     0     0
54     0     0     0     0
55     0     0     0     0
56     0     0     0     0
57     0     0     0     0
58     0     0     0     0
59     0     0     0     0
60     0     0     0     0
61     0     0     0     0
62     0     0     0     0
63     0     0     0     0
64     0     0     0     0
65     0     0     0     0
66     0     0     0     0
67     0     0     0     0
68     0     0     0     0
69     0     0     0     0
70     0     0     0     0
71     0     0     0     0
72     0     0     0     0
73     0     0     0     0
74     0     0     0     0
75     0     0     0     0
76     0     0     0     0
77     0     0     0     0
78     0     0     0     0
79     0     0     0     0
80     0     0     0     0
> 

So in our simulation study, what we actually want is:

> lags_use_multi <- list(R_t = 0:3,
                       S_t = 1:7)

Double check with Carrie when we chat, but this is making more sense now.

carrieholt commented 11 months ago

Just wanted to mention that the lags will be different for R_t and R_prime_t. If we are using R_t=0 in the lags for predicting R_t=t+1, then if we switch back to predicting R_prime_t (recruitment aligned by brood year) ever, S_t=t+1 would be used to predict R_prime_t=t+1 (i.e., in the same brood year). I think this is confusing for me as the Larkin model typically uses notation where the prediction of recruitment is for brood year at t=0, with the independent variable, spawners at t=0, t-1, t-2, and t-3.