Skip to content

incorrect bias correction if bias.correct.control$split is not sorted #11

@brianstock-NOAA

Description

@brianstock-NOAA

Sorry for not taking the time to give a reproducible example... but hopefully this is described well enough for now.

If I specify bias.correct.control$split using Which as in TMBhelper::fit_tmb (lines 130-148), I get incorrect answers if Which is not sorted.

  1. Bias correct all (default):
# mod is optimized TMB model
sdreps <- list()
sdreps[[1]] = summary(TMB::sdreport(mod, bias.correct=TRUE))
  1. Bias correct only a subset of derived quantities, one at a time (very slow, ~4x the above run time)
to.correct <- c("log_SSB", "log_SSB_FXSPR", "log_FAA_tot", "log_FXSPR", "log_pred_catch")
Which = as.vector(unlist(mod$env$ADreportIndex()[ to.correct ] ))
sdreps[[2]] = summary(TMB::sdreport(mod, bias.correct=TRUE, bias.correct.control=list(sd=FALSE, split=Which)))
  1. Bias correct same subset, but hoping to speed it up by using one chunk, split = list(Which)
sdreps[[3]] = summary(TMB::sdreport(mod, bias.correct=TRUE, bias.correct.control=list(sd=FALSE, split=list(Which))))
  1. Bias correct same subset, but sort the indices, split = list(sort(Which))
sdreps[[4]] = summary(TMB::sdreport(mod, bias.correct=TRUE, bias.correct.control=list(sd=FALSE, split=list(sort(Which)))))

Look at bias corrected SSB:

ssb.bc <- lapply(sdreps, function(x) {
	ind.ssb <- which(rownames(x) == "log_SSB")
	bc <- exp(x[ind.ssb,3])
	return(bc)
})
ssb.bc

I find that 1, 2, and 4 give the same bias corrected SSB, but 3 is different (wrong). 2 is very slow, but robust to issues where there are a few NAs (won't return all NAs if there is one NaN). 4 is fast and correct, and should probably be the default in TMBhelper::fit_tmb if BS.control[["vars_to_correct"]] is used. Or maybe TMB::sdreport should check for this?

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions