Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add output for log_pred and log_pred_sd for total biomass and total cpue #16

Open
JaneSullivan-NOAA opened this issue Oct 5, 2023 · 3 comments

Comments

@JaneSullivan-NOAA
Copy link
Collaborator

JaneSullivan-NOAA commented Oct 5, 2023

And maybe the CV on the arithmetic scale for users who want to use total biomass as inputs to another model like SS3?

total_predicted_biomass <- tidyr::expand_grid(model_name = rema_model$input$model_name,

Notes:
Arithmetic-scale CV to log-scale standard deviation (sigma):
sigma=sqrt(log(CV^2+1))
Log-scale standard deviation (sigma) to arithmetic-scale CV:
CV=sqrt(exp(sigma^2)-1)

Example (will need to strip out m1$ in source code):

m1 <- fit_rema(input)

tidyr::expand_grid(model_name = m1$input$model_name,
                   variable = 'tot_biomass_pred',
                   year = m1$input$data$model_yrs) %>%
  dplyr::mutate(log_pred = m1$sdrep$value[names(m1$sdrep$value) == 'log_tot_biomass_pred'],
                sd_log_pred = m1$sdrep$sd[which(names(m1$sdrep$value) == 'log_tot_biomass_pred')],
                log_pred_lci = log_pred - qnorm(0.975) * sd_log_pred,
                log_pred_uci = log_pred + qnorm(0.975) * sd_log_pred,
                pred = exp(log_pred),
                pred_cv = sqrt(exp(sd_log_pred^2) - 1),
                # assumed 95% confidence interval when you use a mean and CV on
                # the arithmetic scale and assume it's normal on the log-scale
                pred_lci = exp(log_pred - qnorm(0.975) * sd_log_pred),
                pred_uci = exp(log_pred + qnorm(0.975) * sd_log_pred)) 

Alternate method (apparently there's no statistical justification to do it one way or the other):

out <- tidy_rema(m1)
out$biomass_by_strata %>% 
  mutate(var_pred = (pred * (sqrt(exp(sd_log_pred^2)-1)))^2) %>% 
  group_by(year) %>% 
  summarise(var_pred = sum(var_pred),
            pred = sum(pred)) %>% 
  mutate(sd_pred = sqrt(var_pred),
         cv_pred = sd_pred/pred,
         sd_log_pred = sqrt(log(cv_pred^2+1)))

@careymcgilliard let me know if you want additional or different output

@Cole-Monnahan-NOAA
Copy link

@JaneSullivan-NOAA I suggest correcting the formula for sigma above which needs the +1 in the log statement

@JaneSullivan-NOAA
Copy link
Collaborator Author

@Cole-Monnahan-NOAA thanks, I thought I did that already...

@Cole-Monnahan-NOAA
Copy link

Ok fixed now, maybe just need to refresh?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants