Modular example of TMB functions

Summary

The goal of this example is to demonstrate how to incorporate key TMB functionalities into a modular C++ framework. The ar1xar1 .R and .cpp code from TMB examples were modified to demonstrate the following TMB functions in a modular framework:

• SIMULATE
• SEPARABLE
• DATA_VECTOR_INDICATOR (for model validation)

File Structure Overview

The ar1xar1 modular example uses the following file structure:

1. Common.hpp

• Includes all header statements
  • #include <TMB.hpp>: includes TMB core files
  • using namespace density: includes functions from TMB’s density.hpp (AR1, SEPARABLE)
• Defines model_traits structure

2. model.hpp: Sets up ar1xar1 C++ templated class

• Line 15 public: defines input data and parameters using the model_traits struct
• Lines 31-33: set up C++ singleton class. This singleton class will allow the functions within to be called within the TMB objective_function
• Lines 36-45: model is defined within the evaluate function

3. ar1xar1.cpp: TMB interface with model.hpp

• Line 2: include model.hpp
• Line 10: set-up TMB objective_function
• Line 13: create a pointer, inst, to point to the singleton class of ar1xar1
• Lines 15-21: set-up data and parameters using TMB macros
• Lines 23-27: Define public members of ar1xar1 class
• Line 29: Call the evaluate function
• Lines 31-37: Simulate data and random effects using the TMB SIMULATE function
• Lines 39-42: Report out derived values and their standard errors

4. ar1xar1.R: I/O R Interface

• Line 4: compiles ar1xar1.cpp
• Lines 6-31: Simulate data
• Lines 37-52: Fit model and generate reports
• Lines 54-58: Simulate data from model and compare to original data
• Lines 60-68: Model validation
• Lines 70-100: Repeat steps for mis-specified model

Incorporating TMB functions into modular framework

1. SIMULATE

The SIMULATE function is run in the TMB interface, ar1xar1.cpp. Simulations are implemented as they would be for a standalone TMB model. For detailed examples, see the comprehensive TMB documentation. If simulations require derived mean values, functions within the singleton class can be called within the SIMULATE function to calculate the mean. See logisticGrowth.cpp for an example. Note that simulation blocks are not commutative, therefore care must be taken to place random effects and data in their proper order.

Data can be simulated from the model from the R interface:

## Not run
sim <- obj$simulate()
plot(y, sim$y)

2. SEPARABLE

Functions within density.hpp, including AR1 and SEPARABLE, are accessible to the ar1xar1 class after adding the following to Common.hpp:

## Not run
using namespace density;

The SEPARABLE fuction requires the random effect vector be of a PARAMETER_ARRAY type. The following code is used to declare the parameter type for the random effect vector, eta, within the model_traits structure of Common.hpp:

## Not run 
tmbutils::array<Type> eta;

The SEPARABLE function is run within the evaluate function of model.hpp.

3. DATA_VECTOR_INDICATOR

TMB uses this MACRO as part of its model validation capabilities for random effects models. For an overview, see the comprehensive TMB documentation.

The DATA_VECTOR_INDICATOR, keep, is defined in the ar1xar1.cpp file as it would in a standalone TMB script. This vector is then used as the definition of the keep vector defined in the model_traits structure of common.hpp and linked through the inst pointer.

//model.hpp

public:
  data_indicator<tmbutils::vector<Type> , Type> keep;

//ar1xar1.cpp

DATA_VECTOR_INDICATOR(keep,y);
inst->keep = keep;

The keep vector can then be used for model validation within the evaluate function of model.hpp:

for(int i=0; i < y.size(); i++){
  nll -= keep[i] * dpois(y[i], exp(eta[i]), true);
  Type cdf = squeeze( ppois(y[i], exp(eta[i])) );
  nll -= keep.cdf_lower[i] * log( cdf );       // NaN protected
  nll -= keep.cdf_upper[i] * log( 1.0 - cdf ); // NaN protected
}

The ar1xar1.R script provides an example of model validation using the original correctly specified model

And using a mis-specified model: