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.

Sign up for GitHub

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

make missing data work with -999 #540

Closed
wants to merge 1 commit into from
Closed

make missing data work with -999 #540

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions inst/include/common/data_object.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ struct DataObject : public fims_model_object::FIMSObject<Type> {
size_t jmax; /**< 2nd dimension of data object>*/
size_t kmax; /**< 3rd dimension of data object>*/
size_t lmax; /**< 4th dimension of data object>*/
Type na_value = -999; /**< specifying the NA value >*/

/**
* Constructs a one-dimensional data object.
Expand Down
30 changes: 20 additions & 10 deletions inst/include/population_dynamics/fleet/fleet.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,11 +168,17 @@ struct Fleet : public fims_model_object::FIMSObject<Type> {
observed_acomp.resize(this->nages);
expected_acomp.resize(this->nages);
Type sum = 0.0;
bool containsNA = false; /**< skips the entire year if any values are NA */
for (size_t a = 0; a < this->nages; a++) {
size_t i_age_year = y * this->nages + a;
sum += this->catch_numbers_at_age[i_age_year];
if(this->observed_agecomp_data->at(y,a) != this->observed_agecomp_data->na_value){
size_t i_age_year = y * this->nages + a;
sum += this->catch_numbers_at_age[i_age_year];
} else{
containsNA = true; /**< sets to true if any values are NA >*/
break;
}
}

if(!containsNA){
for (size_t a = 0; a < this->nages; a++) {
size_t i_age_year = y * this->nages + a;
expected_acomp[a] = this->catch_numbers_at_age[i_age_year] /
Expand All @@ -184,9 +190,10 @@ struct Fleet : public fims_model_object::FIMSObject<Type> {
<< "has expected: " << expected_acomp[a]
<< " and observed: " << observed_acomp[a] << std::endl;
}
dmultinom.x = observed_acomp;
dmultinom.p = expected_acomp;
nll -= dmultinom.evaluate(true);
dmultinom.x = observed_acomp;
dmultinom.p = expected_acomp;
nll -= dmultinom.evaluate(true);
}
}
}
FLEET_LOG << "Age comp negative log-likelihood for fleet," << this->id
Expand All @@ -201,15 +208,18 @@ struct Fleet : public fims_model_object::FIMSObject<Type> {
#ifdef TMB_MODEL
fims_distributions::Dnorm<Type> dnorm;
for (size_t i = 0; i < this->expected_index.size(); i++) {
dnorm.x = fims_math::log(this->observed_index_data->at(i));
dnorm.mean = fims_math::log(this->expected_index[i]);
dnorm.sd = fims_math::exp(this->log_obs_error[i]);
nll -= dnorm.evaluate(true);
if(this->observed_index_data->at(i) != this->observed_index_data->na_value){
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like if we add more data types besides an index and age comps, more if() statements would need to be added to each of the data types' negative log likelihood assignment statements to allow for NAs to work, correct?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

correct! When coding we also discussed if we need to handle a case where there are NA values in some ages but not others despite the existence of some comp samples, but for now, you can only specify a full missing year. Are you thinking we should document this or add to an issue?

dnorm.x = fims_math::log(this->observed_index_data->at(i));
dnorm.mean = fims_math::log(this->expected_index[i]);
dnorm.sd = fims_math::exp(this->log_obs_error[i]);
nll -= dnorm.evaluate(true);
}

FLEET_LOG << "observed index data: " << i << " is "
<< this->observed_index_data->at(i)
<< " and expected is: " << this->expected_index[i] << std::endl;
FLEET_LOG << " log obs error is: " << this->log_obs_error[i] << std::endl;

}
FLEET_LOG << " sd is: " << dnorm.sd << std::endl;
FLEET_LOG << " index nll: " << nll << std::endl;
Expand Down
8 changes: 8 additions & 0 deletions tests/testthat/test-fims-estimation.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -633,8 +633,16 @@ test_that("run FIMS in a for loop", {
catch <- em_input$L.obs$fleet1
fishing_fleet_index <- new(fims$Index, length(catch))
fishing_fleet_index$index_data <- catch
testindex <- 2
na_value <- -999
if(i==4){
fishing_fleet_index$index_data[testindex] <- na_value
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is how R users would use it currently, correct? Assigning -999 to anything that should be an NA?

Maybe this should be added to the PR description at least for some minimal documentation on how this works.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point - I added something to the fims-demo vignette about how to specify missing data!

}
fishing_fleet_age_comp <- new(fims$AgeComp, length(catch), om_input$nages)
fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1
if(i==5){
fishing_fleet_age_comp$age_comp_data[testindex] <- na_value
}

survey_index <- em_input$surveyB.obs$survey1
survey_fleet_index <- new(fims$Index, length(survey_index))
Expand Down
Loading