gadget-framework / gadget3

TMB-based gadget implemtation
GNU General Public License v2.0
8 stars 6 forks source link

error with age-based index #142

Closed vbartolino closed 4 months ago

vbartolino commented 4 months ago

I'm having problem to build a model including an age-based index. This is useful in a number of situations where an index has to be restricted to selected ages, such as a recruitment index, when certain ages are known to be unavailable to the survey or when indices are only available by age as it is the case for many shared stocks via the ICES working groups.

This was working with g2

si.her.aco <- mfdb_sample_count(mdb, c('age'), list(
    area          = mfdb_group('1'="SD30"),
    timestep      = mfdb_timestep_quarterly,
    year          = 2007:2020,
    species       = defaults$species,
    age           = mfdb_group('age1'=1),
    data_source   = 'acoustic_bias_age_her_sd30'))[[1]]

but returns an error in g3 when I build the tmb model

> tmb_model <- g3_to_tmb(actions)
Error in convert_subset(head(cols, -1)) : 
  Missing values must be at start of subset, can't restructure array: adist_surveyindices_log_si.her.aco_model__num[, adist_surveyindices_log_si.her.aco_model__agegroup_idx,     , adist_surveyindices_log_si.her.aco_model__area_idx]

I've tried several alternatives without success including si.her.aco <- mfdb_sample_count(mdb, c('age','length'), list(...

Any suggestion is welcomed, thanks

lentinj commented 4 months ago

The core of the problem is that TMB doesn't know how to do x[, i, , j] (note the gaps between indices). This is arguably for your own good, since chopping up the array into bits like that isn't very efficient, but means we have to be careful with dimension ordering.

What's the definition of the action causing the error? You can put subsets of actions or single actions into g3_to_tmb, e.g.g3_to_tmb(list(g3a_likelihood_distribution(...))).

What does dimnames(g3_distribution_preview(si.her.aco)) look like?

vbartolino commented 4 months ago

What's the definition of the action causing the error?

Not sure I fully understand the question. As I can see the error is caused by the likelihood action and in the specific by the following likelihood component which is associated to the age-based index

  g3l_abundancedistribution(
    'si.her.aco',
    si.her.aco,
    fleets = list(),
    stocks = stocks,
    g3l_distribution_surveyindices_log(beta = 1), 
    area_group = mfdb::mfdb_group('1' = 1),
    nll_breakdown = nll_breakdown,
    report = lik_report)

What does dimnames(g3_distribution_preview(si.her.aco)) look like?

> dimnames(g3_distribution_preview(si.her.aco))
$length
[1] "0:Inf"

$age
[1] "1:1"

$time
 [1] "2007-04" "2008-04" "2009-04" "2010-04" "2011-04" "2012-04" "2013-04"
 [8] "2014-04" "2015-04" "2016-04" "2017-04" "2018-04" "2019-04" "2020-04"

$area
[1] "1"
lentinj commented 4 months ago

Sorry, my comment was a bit cryptic. But thanks for answering regardless!

The error comes from g3l_distribution_surveyindices_log() trying to make a time-series out of the aggregated data, so we can do the linear regression on it. We always put area at the end to get this to work (see the dimnames output). This is a bit of a hack, but it worked up until now.

If $age and $time were also swapped, then the problematic lump of code becomes __num[, , agegroup_idx, area_idx], and we can do what we need.

Regardless, I don't think you're doing anything wrong, we've just not tried this :) I'll have a look tomorrow or Monday to see if we can be more clever about how to order these.