TMB Preliminaries

class: center, middle, inverse, title-slide

# TMB Preliminaries
### Andrea Havron NOAA Fisheries, OST

---

layout: true

.footnote[U.S. Department of Commerce | National Oceanic and Atmospheric Administration | National Marine Fisheries Service]

code.cpp{
  font-size: 14px;
}
code.r{
  font-size: 14px;
}

</style>

---

# TMB Model

.pull-left[
[**linReg.cpp**](https://github.com/NSAWTraining/TMBtraining/blob/main/docs/01-beginner-tmb/lecture_code/linReg.cpp)
```cpp
#include <TMB.hpp>
template <class Type>
 
Type objective_function<Type>::operator()()
{
 DATA_VECTOR(y);
 DATA_MATRIX(X);
 PARAMETER_VECTOR(beta);
 PARAMETER(lnSigma);
 
 Type nll = 0;
 Type sigma = exp(lnSigma);
 int n = y.size();
 vector<Type> mu = X * beta;
 for(int i=0; i<n; i++){
 nll -= dnorm(y(i), mu(i), sigma, true);
 }
 SIMULATE{
 y = rnorm(mu, sigma);
 REPORT(y);
 }
 
 Type sig2 = pow(sigma,2);
 REPORT(sig2);
 ADREPORT(sig2);
 
 return nll;
}
```
]

.pull-right[
**C++ Syntax**
* Lines end in semi-colons
* Everything must be declared
* Type is generic and assigned based on input
* Indexing starts at `$0$`!
* `$x -= 1: x = x-1$`
* Matrix multiplication uses `$*$`
* Math operators similar to R `$(+ - /*)$` 
* Use `$pow(x,p)$` for `$x^p$`
* if statements cannot be based on a parameter
]

---
# Data: Importing data from R 
* Pass data to TMB with these 'macros'
* Note: do not specify the object dimension

|TMB Syntax |C++ Type |R Type |
|:---------------|:----------------------|:----------|
|DATA_VECTOR(x) |tmbutils::vector<Type> |vector |
|DATA_MATRIX(x) |tmbutils::matrix<Type> |matrix |
|DATA_SCALAR(x) |Type |numeric(1) |
|DATA_INTEGER(x) |int |integer(1) |
|DATA_FACTOR(x) |Eigen::vector<int> |factor |
|DATA_ARRAY(x) |tmbutils::array<Type> |array |
---
# Data: Importing data from R

.pull-left-narrow[
TMB code
```cpp
DATA_VECTOR(y); 
DATA_MATRIX(X); 
DATA_INTEGER(i); 
DATA_FACTOR(ngroup);

```
]

.pull-right-wide[
R script

```r
Data <- list(
 y = c(30.2, 45.3, 12.1),
 X = matrix(0,3,3),
 i = 11,
 ngroup = c(1,1,2,2,2)
)
str(Data)
```

```
## List of 4
##  $ y     : num [1:3] 30.2 45.3 12.1
##  $ X     : num [1:3, 1:3] 0 0 0 0 0 0 0 0 0
##  $ i     : num 11
##  $ ngroup: num [1:5] 1 1 2 2 2
```
]

---
# Declaring model parameters 
* No PARAMETER_INTEGER
* Again, do not specify the object dimension

|TMB Syntax |C++ Type |R Type |
|:-------------------|:----------------------|:----------|
|PARAMETER_VECTOR(x) |tmbutils::vector<Type> |vector |
|PARAMETER_MATRIX(x) |tmbutils::matrix<Type> |matrix |
|PARAMETER_ARRAY(x) |tmbutils::array<Type> |array |
|PARAMETER(x) |Type |numeric(1) |

---
# Declaring model parameters

.pull-left-narrow[
TMB code
```cpp
PARAMETER_VECTOR(beta); 
PARAMETER(lnSigma); 
PARAMETER_MATRIX(u); 
```
]

.pull-right-wide[
R script

```r
Pars <- list(
 beta = c(0,0),
 lnSigma = 0,
 u = matrix(0,3,3) 
)

str(Pars) 
```

```
## List of 3
## $ beta : num [1:2] 0 0
## $ lnSigma: num 0
## $ u : num [1:3, 1:3] 0 0 0 0 0 0 0 0 0
```
]
---
# Parameter Transformations
.pull-left[
Transform parameters from real to parameter space
 * Minimizers search for MLE from `$-\infty$` to `$+\infty$`
 * Certain parameters have restricted support:
 * `$\sigma^{2} > 0$`
 ```cpp
 sigma = exp(lnSigma)
 ```
 * `$0 < \pi < 1$` 
 ```cpp
 pi = 1/(1 + exp(-logit_pi)) # "expit"
 ```
]

.pull-right[
![](data:image/png;base64,#00_01_TMB_Preliminaries_files/figure-html/unnamed-chunk-5-1.png)![](data:image/png;base64,#00_01_TMB_Preliminaries_files/figure-html/unnamed-chunk-5-2.png)
]
---
# Calling a TMB model from R
 
**Running Optimization**

```r
obj <- MakeADFun(data = Data, parameters = Pars, DLL = 'linReg')
opt <- nlminb(obj$par, obj$fn, obj$gr)
```

**Reporting Results**

.pull-left[
C++
```cpp
  Type sig2 = pow(sigma,2);
  REPORT(sig2);
  ADREPORT(sig2);
```
]
.pull-right[
R

```r
report <- obj$report()
sdr <- sdreport(sdr)
summary(sdr, "fixed")
summary(sdr, "report")
```
]
---
# Parameter mapping

#### TMB allows users to collect and fix parameters using the map argument in MakeADFun()
.large[.p[
* The parameter map is a list of parameters that are fixed or collected in a model run
* Names and dimensions of parameters in the map list must match those of the parameter list
* The map list is structured as a list of factors
* Parameters with factor(NA) are fixed
* Parameters with equal factor levels are collected to a common value
]]

---
# Likelihood Ratio Tests

.pull-left[
`$H_{0}$`: Null Hypothesis 
.p[
- `$\theta \in \Theta_{0}$`
- parameter `$\theta$` is in a subset of a given parameter space, `$\Theta_{0}$`
- restricted model
]]
.pull-right[
`$H_{1}$`: Alternative Hypothesis

.p[
- `$\theta \in \Theta$`
- parameter `$\theta$` is in the complete parameter space, `$\Theta$`
- unrestricted, full model
]]

**Likelihood Ratio Test**:
`$$LR = -2[\ell(\theta_{0}) - \ell(\hat{\theta})]$$`

`$$LR \sim \chi^{2}(df_{\chi^{2}} = df_{\Theta}-df_{\Theta_{0}})$$`

---
# Parameter mapping with LR test
 
TMB example: [**lr_test.cpp**](https://github.com/NSAWTraining/TMBtraining/blob/main/docs/01-beginner-tmb/lecture_code/lr_test.cpp)

```cpp
/ Illustrate map feature of TMB to perform likelihood ratio tests on a ragged array dataset.
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
 DATA_VECTOR(obs);
 DATA_FACTOR(group);
 
 PARAMETER_VECTOR(mu);
 PARAMETER_VECTOR(sd);
 
 Type res=0;
 
 for(int i=0;i<obs.size();i++){
 res -= dnorm(obs[i],mu[group[i]],sd[group[i]],true);
 }
 
 return res;
}

```
---
# Simulation in TMB

* Standard generator functions to simulate data within the TMB model:

```cpp
rnorm()           rpois()           runif()
rbinom()          rgamma()          rexp()          
rbeta()           rf()              rlogis()  
rt()              rweibull()        rcompois()
rtweedie()        rnbinom()         rnbinom2()
```

* Simulation blocks are used to call simulations from R

.pull-left[
**Binomial Example** 
C++ code
```cpp
for(int i = 0; i < n_obs; i++){
 eta(i) = beta(0) + beta(1)*x(i); 
 p(i) = 1/(1 + exp(-eta(i)));
 nll(i) = -dbinom(Y(i), N(i), p(i), true); 
 SIMULATE{
 Y(i) = rbinom(N(i), p(i));
 }
}
SIMULATE{
 REPORT(Y);
}
```
]

.pull-right[
 
R code

```r
set.seed(1) ## optional

obj$simulate(complete = TRUE)
```
*Note*
* optimization not necessary for simulation
]

---
#Simulation Study

.pull-left[
[**linReg.cpp**](https://github.com/NSAWTraining/TMBtraining/blob/main/docs/01-beginner-tmb/lecture_code/linReg.cpp)
```cpp
 for(int i=0; i<n; i++){
 nll -= dnorm(y(i), mu(i), sigma, true);
 }
 SIMULATE{
 y = rnorm(mu, sigma);
 REPORT(y);
 }
```
Set-up Model in R

```r
library(TMB)
dyn.load(dynlib('lecture_code/linReg'))
sig <- 2
beta <- c(0.5,2)

Data <- list(y = rep(0,50), 
 X = cbind(rep(1,50), 1:50))
Pars <- list(beta = c(0.5,2),lnSigma = log(2))
obj <- MakeADFun(data = Data, 
 parameters = Pars, 
 DLL = 'linReg')
obj$par
```

```
##      beta      beta   lnSigma 
## 0.5000000 2.0000000 0.6931472
```
]
.pull-right[
Run Simulation

```r
set.seed(1)
sim <- replicate(500, {
 simdata <- obj$simulate(par=obj$par, complete=TRUE)
 obj2 <- MakeADFun(simdata, Pars, 
 DLL="linReg", silent=TRUE)
 nlminb(obj2$par, obj2$fn, obj2$gr)$par
})
obj$par
```

```
##      beta      beta   lnSigma 
## 0.5000000 2.0000000 0.6931472
```

```r
apply(t(sim), 2, mean)
```

```
## beta beta lnSigma 
## 0.5143324 1.9995115 0.6649257
```
]
---
#Random Effects Models

### The Hierarchical Model:

`$$\Large \int_{\mathbb{R}}f(y;u,\theta)f(u;\theta)du$$`
---
# The Laplace approximation

Changes the problem from integration to optimization 
.pull-left[
`$$L(\theta) = \int_{\mathbb{R}}f(y;u,\theta)f(u;\theta)du$$` 
`$$L^{*}(\theta) = \sqrt{2\pi}^{n}det(\mathbb{H})^{-1/2}f(y,\hat{u}, \theta)$$`
]
.pull-right[
<img src="data:image/png;base64,#static/laplace-accuracy.png" width="100%" style="display: block; margin: auto auto auto 0;" />
 
Figure from [Albertsen, C. M. (2018), 2.3.1](https://backend.orbit.dtu.dk/ws/portalfiles/portal/157133664/Publishers_version.pdf)
]
---
# Running from R

.pull-left[
[**first_re.cpp**](https://github.com/NSAWTraining/TMBtraining/blob/main/docs/01-beginner-tmb/lecture_code/first_re.cpp)
```cpp
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
 DATA_VECTOR(y); 
 DATA_SPARSE_MATRIX(B); 
 DATA_SPARSE_MATRIX(X); 
 PARAMETER_VECTOR(u); 
 PARAMETER_VECTOR(beta); 
 PARAMETER(lnSDu); 
 PARAMETER(lnSDy);

Type nll = 0;
  Type sdu = exp(lnSDu);
  nll1 -= sum(dnorm(u, Type(0), sdu, true));

vector<Type> mu = X * beta + B * u;
 Type sdy = exp(lnSDy);
 for(int i=0; i<y.size(); i++){
 nll2 -= dnorm(y(i), mu(i), sdy, true);
 }

return nll;
}
```
]

.pull-right[

```r
Dat <- list(y=y, B=B, X=X)
Par <- list(u=u*0, beta=beta*0, 
 lnSDu=1, lnSDy=1)
obj <- MakeADFun(Dat, Par,
 random="u", 
 DLL="first_re"
 )
opt <- nlminb(obj$par, obj$fn, obj$gr)
```
]