Closed Lewis-Barnett-NOAA closed 5 months ago
More details on how I envision simulating from the linear model of sea ice on cold pool. You could make the gradient of fish density dependent on cold pool extent as a covariate (when cold pool is low, less gradient in density; when cold pool is high, more gradient in density). Then....
load the fitted model object by providing the path to where you save the .RDS file within your directory:
m <- readRDS("~/Data/ice_coldpool_lm.RDS")
set the seed to an arbitrary number so that you can replicate your results from a given simluation instance if you start with that same seed:
set.seed(7)
sims <- simulate(m, nsim = 10)
This will give you, for example, 10 replicates of simulated cold pool extents (but you can specify as many as you want, default is 1), which correspond in order to each of the values of march sea ice that went into the initial model fit, which you can obtain from this object:
m$model$march_sea_ice
For example, you can plot the observed sea ice against simulated cold pool for any given simulation as:
plot(m$model$march_sea_ice, sims$sim_1)
So, now you can combine the observed sea ice values with as many simulated cold pool extents resulting from each sea ice value:
cbind(m_lte2_prop$model$march_sea_ice, sims)
In your simulation routine, for each iteration you would then draw a random value of sea ice from the first column to use as your indicator of whether to sample more or less in the strata that typically has fewer observations. You would also draw the corresponding prediction of cold pool extent from one of the columns of simulated responses. If cold pool is used as a covariate driving the density gradient, you can then plug in the simulated cold pool extent value into your model for fish density and simulate survey samples from that surface with observation error. Or, you could just arbitrarily assign different gradients of density (different spatial models) to be your representation of the true density distribution given different levels of cold pool, and sample from those in the same way.
Here the procedure would be the same as in stage 1, but instead of the canned simulate() function, you could use your own function to perform the simulation from the linear model of sea ice and cold pool. This code replicates what simulate() does:
for(i in 1:nsim) {
tmp <- predict(m) + rnorm(length(predict(m)), 0, summary(m)$sigma)
res <- if (i==1) tmp else cbind(res, tmp)
}
head(res)
So, now you can increase or decrease the noisiness of the predictions (or inversely the predictive skill of our indicator) by increasing or decreasing the variance from the fitted model, for example, by changing summary(m)$sigma
to 0.1*summary(m)$sigma
in the above code, you would get simulated cold pool values with only 10% of the estimated standard deviation from the originally fitted model.
The two stages above is the general approach that I recommend for now, but I'll outline one more option for future reference.
You can modify the simulate() function to give simulated cold pool values for any value of sea ice we could observe (lets say 0 to 1 by 0.05). This way we can explore the full range of possible values for sea ice and cold pool. Here is a link to two ways of doing this: https://stackoverflow.com/questions/14967813/is-there-a-function-or-package-which-will-simulate-predictions-for-an-object-ret
All in this repository, see the .R script for details about the sea ice cold pool regression: https://github.com/danielvilasgonzalez/Bering_redesign/tree/main/Data/Sea_ice_data
Really, all we need for this project is ice_coldpool_lm.RDS which is the fitted linear model, which we can simulate from to serve as an indicator. We can also use this model of an indicator to explore a range of predictive skill of an indicator, by changing the amount of observation error (basically inflating the variance, or decreasing the R^2, and vice versa).