sdmTMB Reduction Simulation

[AngeliaMiller]

Hi @gavinfay,

I tried looking into the results for the productivity scenarios based on the treatment effects. Specifically, there was a question about the reduction scenario and its smaller relative percent differences when lower productive wind areas are precluded from the survey index calculation. To derive the plot included in the ICES presentation, I forced the wind area coefficient to be 0.09-log(2) for the reduction scenario.

I tried changing the wind energy coefficient to 0.09-log(1/2) and simulating only 25 replicates for time to receive the resulting plot, but I think this would be the same as 0.09+log(2) if I understand my log rules correctly.

Next I tried 0.09+log(1/2) and simulated 25 replicates and received the resulting plot:

The latter plot seems a little better, but still not quite sure what I'm looking for to properly diagnose the issue. The distribution of estimates of abundance under these two trials are also not seeing an increase when you remove the lower productive areas from the sample frame.

Any feedback you could provide would be greatly appreciated.

Thanks!

[gavinfay]

@AngeliaMiller We talked about not trying to diagnose things so far from where the problems are likely arising. This plot is pretty far from the data, it is not only a summary statistic over years from several simulations, but it is also a comparison among scenarios.

You don't need to revisit regression coefficient values, you already know these are correct, i.e. e^(reduced) = 0.5 & e^(enhanced) = 2.

Thinking about what we expect to see in this plot: The plot is intending to show the distribution over simulations of the differences between abundance estimates when tows from wind areas are included or not in the calculation. We expect these differences to be smallest (closer to 0) for the scenario with the smallest treatment (the smallest difference between the expectation for catch rates in and outside the wind areas), which is the baseline scenario. While the math underlying it is not correct, because both scenarios are the same enhanced scenario, the 2nd of the 3 plots in your post above is more in line with what we would think we would see. That doesn't mean that we will see that result (if we knew exactly what would happen there would be no need to do the research), but when something is very contrary to expected, we need to figure out why.

As mentioned in a previous email, it could be that your mean is low enough (particularly in the reduced scenario) that under the tweedie you are in the space of a poisson where most of the probability mass is piled up so there isn't much difference between the probability distributions even when the means are very different, whereas under the baseline you do at least have some observations that results in an extended upper tail, thus meaning the overall abundance estimate is sensitive to inclusion or not of these values (though I find it slightly surprising the overall abundance estimate is so sensitive to this). Perhaps you would get the expected result if the scale was different and you had a higher mean. This sort of thing you could get a feel for very quickly by testing in R with random samples from distributions with various sets of parameters (ie no need for the survey simulation code).

An aside (sort of): What is being plotted on the y axis here? It seems unlikely that it is the distribution of percent relative differences as indicated in the axis label, maybe the distribution of absolute percent relative differences?

So I might proceed by (might not have to do all of these depending on whether you figure out if/where the issue lies):

1. checking the summary statistic being presented in the graphs here.
2. doing the check mentioned above by generating distributions from a tweedie (or even a gamma or poisson, depending on what value your tweedie power parameter is), with and without your treatments (reduction coefficient), comparing the differences and then do the same with alternatives that change the scale of the observations by increasing the underlying mean. This should give you a feel for whether (if all tows are considered equal) you should be seeing what you think.
3. Step through the analysis just by considering 1 simulation of 1 year and
4. Verify that the catch rates from the simulated survey tows occurring in wind areas have the intended treatment effect (compare the distribution of tow catch rates in and outside areas, controlling for covariates (i.e. maybe use tows from the same strata as these should have similar depths). You might want to run a few replicates (either via simulations or years) here so you can get sufficient sample size to look at (as the number of tows in wind areas in a given survey is probably low).
5. Check your abundance mean and variance calculations. e.g. Pick a stratum that contains wind areas, and check that the stratum mean and variance are being calculated right and that they are sensitive to the preclusion scenario in the right way. e.g. copy the data from these twos into Excel and do the calculations by hand, and compare to your R function.
6. Check your overall annual abundance mean and variance calculations - are these being done correctly over strata? Do this (and 4) with dummy data with just a few (say 3) strata and observations within them, with the 1 stratum containing wind areas and the others not. i.e. create a simple simulation experiment of the process that you can test your functions with rather than getting bogged down with the complexities of the full data set and survey strata. As stated in the past, a test toy data set version for this will be useful for testing many other parts of this project, and a more comprehensive simulation study that is generic but conditioned on the use case would form a great basis for a paper/dissertation chapter, very much in the spirit of Debie Duarte's initial PhD chapter with her simulation experiment looking at observer effect.

[AngeliaMiller]

Hi @gavinfay,

I worked through steps 2-4 delineated above for both summer flounder and Atlantic mackerel and generated some respective plots. Here I am attaching the plots generated based on the fall summer flounder model fit; the respective code can be found here. Plots for the spring summer flounder model and spring Atlantic mackerel model can be found here beginning on slide 53 for summer flounder and slide 103 for Atlantic mackerel. Let me know how I can adjust the plots in anyway to make them more informative.

Thanks!

---

2. Generate distributions from a tweedie with and without treatments (e.g increase the mean by log(2) and decrease the mean by log(2)). -- the phi and power parameters estimated by the fall summer flounder model were used to generate a distribution of 100 numbers around a mean of 1, 10, and 100 -- each mean was then increased or decreased by log(2) to reflect the intended treatment for our productivity scenarios

3. Step through the analysis just by considering 1 simulation of 1 year -- plot reflects the 2016 annual abundance index for summer flounder given the survey effort scenario (status quo or preclusion) and the wind-driven productivity scenario (increased or decreased catch rates in wind energy areas).

5. Verify that the catch rates from the simulated survey tows occurring in wind areas have the intended treatment effect -- 5 replicates were simulated for all years at the same tow locations as those to fit the model -- each replicate was filtered for a given wind-impacted strata, and where catch rates were greater than 0. -- the plot only shows the distribution of wind energy area catch rates within this stratum to verify the intended treatment effect

[gavinfay]

Ok, so what are your conclusions? Noting that you haven't yet compared the differences here. I'm not sure the simulated random numbers are really including your treatment effect as intended - e.g. for the case when the mean is 100, the 'increase' scenario clearly doesn't have a mean of 200.

[AngeliaMiller]

My conclusions are that our treatments are having the intended effect - reducing the estimated wind coefficient will simulate decreased catch rates in wind energy areas of summer flounder, and increasing the estimated wind coefficient will simulate increased catch rates in wind energy areas of summer flounder.

Additionally, based on the simulated random numbers, a lower mean with the fitted tweedie distribution results in a piled probability mass, and given a further reduction of an already small mean (eg. the reduction treatment), there would be no difference in the probability distribution. Thus, in the summer flounder case, the abundance index will be more sensitive under the reduced scenario, than it would be under a baseline effect. (I'm not sure I summarized that correctly)

With regards to your last point, I simulated the random numbers such that the mean = 100 + log(2), so if we wanted a mean of 200, I guess I should have simulated it such that the mean = 100*2 ?

[gavinfay]

With regards to your last point, I simulated the random numbers such that the mean = 100 + log(2), so if we wanted a mean of 200, I guess I should have simulated it such that the mean = 100*2 ?

You're mixing up response space and link space. The treatment in your experiments is to double (x2) catch rate inside wind areas when productivity is increased, and to halve it (x0.5) in the decreased producitvity scenario.

Your statements above are in line with the hypothesis (sort of, you're focused on effect size rather than comparison of effect size to baseline), but you haven't demonstrated these conclusions yet.

For item (2), first implement the correct effect size, and then you want to look at the sampling distribution for the differences between data from the perturbed sites (with your effect) and unperturbed (with no effect), for all three effect scenarios. You can get a sense for the results by just plotting the distributions of the initial random draws but this isn't getting at the issue of concern fully (whether there is an expected reduction in differences compared to baseline when the effect is the reduction scenario when the mean catch rates are low).

For item (3) it would be good to go a step further and calculate the absolute percentage difference. For what you show in the plot, it looks like your initial results (with a smaller difference in the reduction scenario) are wrong - here you show the expected larger effect in the reduction scenario than under baseline. (though granted this is just 1 simulation)

For item (4), you aren't showing that the survey data has the intended treatment effect - you're only showing data from inside the wind areas, not comparing tows inside and outside the wind areas in the same stratum. The comparison is inside/outside wind given a treatment (or baseline), not within wind areas across treatments.

Hope this helps.

[AngeliaMiller]

For item (2) I corrected the effect size. To look at the sampling distribution for the differences between data from the perturbed sites and the unperturbed sites, I interpreted this as drawing 1000 random samples of size 50 from the initial randomly generated 1000 numbers under each effect and found the mean of each sample. The distribution of 1000 means is plotted here. This piece here:

whether there is an expected reduction in differences compared to baseline when the effect is the reduction scenario when the mean catch rates are low

was tripping me up a little bit in terms of how to enact your suggestion.

For item (3), here is the absolute percentage difference for that 1 simulation. Calculated as:

• the difference between the survey effort scenario stratified means for each productivity scenario
• the average of the survey effort scenario stratified means for each productivity scenario
• the difference divided by the average, multiplied by 100.

For item (4), I made a few more plots that compare the catch rates inside and outside wind areas inside the same stratum given a treatment.

[AngeliaMiller]

Hi @gavinfay,

Following up from our meeting on Friday regarding this troubleshoot.

I was able to generate tweedie distributions of random numbers from the fall power and phi parameters around three scales of means (1, 10, and 100), with and without the treatment effect, and compare the sampling distribution differences (e.g. item 2 above)

This exercise helps confirm that the baseline scenario should have the smallest differences, as was expected.

Then, I generated 5 iterations of productivity changes only in a strata impacted by wind for all the years it was surveyed, calculated the mean in a status quo effort and preclusion effort scenario for each year and replicate, found the absolute percent difference between the two means for each year and replicate. However, when I plot the distribution of these differences (5 replicates x 11 years = 55 differences), I'm not seeing the same intended effects, specifically when looking at the enhanced scenario, whose median difference is the same as the baseline, but lower than the reduction scenario. Would this because we are looking at the extreme case here? It seems like even in the extreme case, the enhanced scenario should still have the largest difference, followed by the reduction scenario, and then the baseline scenario with the smallest difference.

[AngeliaMiller]

Hi @gavinfay,

Based on our meeting on 01/25, I was able to simulate changes in one wind-overlapped stratum and year and calculate the absolute percent differences between survey effort indices for both fall and spring summer flounder distributions. I was also able to simulate the null hypothesis to look at when there is no wind effect on the catch rates for summer flounder. I put all the plots on this Google slide deck. I will note that it did take a few combinations of strata and year before I finally saw the expected results. Certain combinations were still giving me results where the differences for the baseline were higher than the reduction scenario. All my code can be found here

Similarly, I was able to conduct all the same checks for Atlantic mackerel except the one stratum and one-year analysis. I have tried about 20 combinations of year and stratum and each time, the differences for the baseline were either greater than or equal to the reduction scenario differences. I'm wondering if this is because Atlantic mackerel is so variable in its observations regardless of whether catch rates occur in wind areas or not, but then that doesn't line up with the expected loss analysis that the losses to catch rates and the abundance index are 15% and 30% respectively when survey effort is reduced. Conversely, though, I am seeing the intended effect in the CVs between productivity scenarios, such that the baseline has the lowest differences compared to the enhanced and reduced scenarios.

[gavinfay]

@AngeliaMiller looking at this quickly it seems you are calculating means and variances for each simulation? The goal of doing multiple replicates was to increase samples size because there are so few tows in wind areas in a given stratum in a given year. I.e. to remove the effect of sampling variability. That way you would get a better sense of the distribution of catch rates from both inside and outside the wind area so you can see if your intended treatment effect is being generated correctly. I am not surprised that your results are variable if you are not increasing n but doing a bunch of simulations with a small number of samples.

[AngeliaMiller]

@gavinfay,

Thank you for the feedback. I've updated the slide deck to include the distribution of catch rates within one year and strata across 500 simulations of our three treatments for both our species and their respective seasons.

Based on the resplots in this slide deck, it does seem like we are generating the expected treatments.

Revisiting your troubleshooting steps above, all that remains are steps 5 and 6, which may be addressed by the issue you recently found regarding the relative strata weight and total survey area application within the stratified mean functions, but something that I intend to look into as a result.

Thanks!

[gavinfay]

@AngeliaMiller I am not seeing the plotted distributions of simulated catch rates that you refer to. Can you please point me to them? Thanks.

[AngeliaMiller]

@gavinfay, sorry about that. They are on slides 7, 15, and 24. I have removed the zeros from the distribution plots as well. Thanks!

[gavinfay]

OK, I see what you are referring to now.

Based on the resplots in this slide deck, it does seem like we are generating the expected treatments.

I don't see from these how you can tell?

Some challenges here:

• you are collapsing (almost) all the observations in each distribution into a single bin, so no real way of comparing
• Something appears wrong with the simulated data. The range of observations in your simulated survey seems off (0 - 60,000 kg per tow; do we see individual survey tows catching 60 tonnes of summer flounder?). Simulated catch rates for this same strata (1650) for summer flounder earlier in this thread were in the range 0-20 kg per tow.
• why remove the zeroes? They are part of the observed distribution, and one of the reasons you had postulated for not being able to detect an effect in your alternative senarios.
• assuming you fix the binning, comparing distributions here is hard because they are offset from eachother due to being overliad, this hard to make direct visual comparison. I would have expected to see something like the below.

[AngeliaMiller]

Hi @gavinfay

Thanks for the feedback. I am looking into the simulated catch rates. I will note that I did find my labels are wrong on the plots, specifically for summer flounder; the simulated catch rates in the slide deck for the fall pertain to catch rates in strata 3200 and year 2018, not 2009 and 1650.

Regardless of that, you are correct that something is wrong with the simulated catch rates. In stratum 3200 the historically observed catch rates range from 0 - 3 kg/tow; so we do not see individual tows with 60 tonnes of summer flounder. Similarly, when I simulate catch rates for fall summer flounder in strata 1650, I am still generating large values of catch rates, with the maximum being 853 kg/tow, which still goes against the simulated distribution plotted above where it only ranged from 0 - 20 kg/tows; and historical observations in this stratum also only range from 0 - 3 kg/tow.

[AngeliaMiller]

Hi @gavinfay

It seems there was an issue with the way I was extracting and calling my beta values of estimated coefficients in the simulation function. I didn't exactly diagnose what the issue was (e.g my extracted year coefficient did match my selected year so perhaps something with my formula), other than when I used the full set of estimated coefficients to simulate at all of the locations used for the model fit, and then filtered for a given year and strata post-simulation, the catch rates were no longer on the orders of 10^5, and instead were in ranges of 0 - 50 kg/tow for summer flounder (which is how I originally generated the plot you referenced above).

I've updated the slide deck with the distributions of simulated catch rates on slides 6, 13, and 21 generated from 100 iterations. I updated the plots to include simulated catch rates of zeroes as well.

[gavinfay]

Ok. Good that you fixed the error. Can you change the bin widths and the x axis scale on the plots so you can view the data? (ie have the data make up most of the plot rather than being squished up on the very left and lots of empty space). Particularly so for spring fluke and mackerel. Fall fluke looks like it is doing what you expect with the treatments.

[AngeliaMiller]

I've truncated the x-axis's on both the spring fluke and mackerel plots here on slides 13 and 21.

[gavinfay]

OK, so this gives you a bit more of a window in to what these distributions look like (though you can do more with bin width). What are your thoughts given these results? (and also in context with the directly simulated values summarized earlier in the thread - you might find looking at the distribution for the tows there helpful in addition to the distributions of the differences)

[AngeliaMiller]

Here are my thoughts:

• We are seeing the intended effects of our treatments on simulated fluke catch rates in wind areas -- there is a larger quantity of zero catch rates (or close to zero) being simulated in the reduced scenario in each season both compared within wind areas only and comparing wind catch rates to the outside catch rates in the reduced scenario. -- if wind areas are associated with higher catch rates in the baseline, then reducing those catch rates would make them less than or similar to the outside areas ---- Then, the mean catch rate of status quo effort in the reduced scenario would be lower than the mean catch rate in the baseline scenario -- removing reduced catch rates in the preclusion scenario would also mean a higher average catch rate than when they are included ---- the precluded average catch rate in the reduced scenario would be (should be) the same as the precluded average catch rate in the baseline scenario -- so there should be a larger difference between indices of survey effort scenarios in the reduced fish density treatment than the baseline density

• Mackerel on the other hand has almost an equal distribution of catch rates in wind areas across all three scenarios if only slightly greater quantities of lower catch rates in the reduced scenario compared to the baseline -- so while following the same thought process above, there might still be a larger difference between indices of survey effort in the reduced fish density treatment than the baseline, but the magnitude of difference between scenarios might be smaller for mackerel when compared to the magnitude of difference between scenarios for fluke.

[gavinfay]

Good. I think your logic is sound for summer flounder here. So now you know that with adequate sample size, you should see a greater difference in abundance estimates when tows are precluded vs not under your two treatment scenarios compared to the baseline with the estimated coefficients. Now to check your 2-stage sample means/variance estimators as discussed above.

For mackerel, I think your conclusions are still hampered a bit by using large bins in your frequency distributions. But, one immediate question is what is driving the larger proportion of zero tows in the wind areas? Is there a negative regression coefficient for the wind area covariate? The plots would suggest so? You might find it helpful to report the proportion of zeroes along with these, given the continuous nature of the data.

[AngeliaMiller]

Thanks Gavin. The Atlantic mackerel model does not have an estimated area effect from the model fit because model diagnostics suggested not to include it in the linear predictor. So currently when simulating data for the baseline scenario, I include Area in the linear predictor and set its coefficient to 0; a similar approach with with fix_omegas.

I'll work on the bin width and the proportions of zeroes for the one year, one stratum, multiple replicates exercise.

[gavinfay]

OK so something else is causing a difference between the wind area and outside. Did you compare the depth effects for the relevant cells in this stratum as discussed previously? I can't imagine they will be too different because the fluke estimates aren't that different. So the other thing could be the spatial and spatiotemporal coefficients for these cells for the given year you are looking at. I would look there.

[AngeliaMiller]

If by "depth effects" you mean the depth at which the tows occurred previously that we are now simulating a response at, then you are correct, the depths are not that much different between the outside area and the wind areas. If this is not what you mean, then I think I need some additional clarification as to what you mean by comparing depth effects.

With regards to the spatial and spatiotemporal coefficients, do you mean the simulated values of the spatiotemporal random fields, and the model expected values of the spatial random effects (e.g the outputs from sdmtmb_simulate()), or the model expected values of both of these terms (e.g. outputs from predict())? Am I looking at how much they differ between outside and inside wind areas similar to the comparison of depth effects in the given year and strata?

[gavinfay]

Your expected value for each grid cell is something like:

log(u{cell, year}) = year + f(depth) + omega{cell} + eta_{cell, year}

If the year effect and depth effect are the same then the only other things that can be changing are the spatial effect (omega) or the spatiotemporal effect (eta).

[gavinfay]

If it is the spatial effect that is causing the difference then this will persist for this stratum over years. If it is the spatiotemporal effect then the difference you see in the catch rates for this stratum could vary depending on the year, meaning it is important to capture that variability when understanding survey performance. The implications are different. Either way, important to understand what is driving the differences/behavior you are seeing (this is generally true).

[AngeliaMiller]

Okay, so if I am understanding correctly:

I looked at the model-predicted values for omega_s and epsilon_st at the historical survey locations for spring mackerel for 2009 and in stratum 1690. Then I looked at these same columns but at the survey locations in 2010 and stratum 1690. Both epsilon_st and omega_s differ between years, which makes sense because the station locations differ between years.

Then I looked at the model-predicted values for epsilon_st and omega_s at the grid cells for 2009 and stratum 1690, and then again for 2010 and 1690. In this case, the omega_s is the same each year (e.g. the predicted values persist in the stratum over years) and epsilon_st is changing each year.

All the plots and simulations conducted for this thread have been simulating observed biomass at the same station locations as the historical data and thereby fixing the omega_s values with the values predicted from the model at these same locations. I have not simulated biomass in any of the grid cells yet, though I do have the code ready to go to include in the loop.