1 Preparing `hhh4` Data for `RTMB`

As example data we use the demo data from the hhh4contacts package. The data consist of norovirus gastroenteritis counts stratified by age group and district in Berlin from 2011–2015 reported to the Robert Koch Institute (RKI), obtained from SurvStat@RKI.

We use these data since they allow us to later extend the models to include contact matrices and age-group stratification. The data are provided as an sts object and can be aggregated either by district only or by district and age group using the function hhh4contacts::noroBE().

To use the data for RTMB models, we extract the observed counts from the sts object together with the neighbourhood matrix and population data.

# Norovirus data stratified by district
noroBE <- noroBE(by = "districts")

# Extract data from sts object as matrices
noro <- observed(noroBE)
neigh <- neighbourhood(noroBE)
pop <- population(noroBE)

n_time <- nrow(noro)

For the use in RTMB models we create a list containing all required data objects.

# Create data list to use in RTMB models
data <- list(
  Y = noro, # response
  n_time = nrow(noro),
  n_unit = ncol(noro),
  pop = pop
)

2 Endemic-only Models

First we are going to build a simple endemic-only model.

2.1 Mathematical Formulation

In a pure endemic model, we assume the disease incidence in unit $i$ at time $t$ follows a Negative Binomial distribution where the mean $\mu_{it}$ depends only on a baseline endemic component and a population offset. Since there is no autoregressive or neighborhood component, the model simplifies to:

\[Y_{it} \sim \text{NegBin}(\mu_{it}, \psi)\]

\[\mu_{it} = e_i \nu\]

\[\log(\nu) = \beta_0\]

Where $e_i$ is the population size (included as an offset), $\beta_0$ is the endemic intercept, and $\psi$ is the overdispersion parameter. The population offset ensures that the expected number of cases scales with the size of the population. This is structurally identical to a standard Negative Binomial Generalized Linear Model (GLM).

2.2 Using `surveillance`

Using hhh4() we need to create the control list first.

## Control for endemic-only model
noroModel <- list(
  # endemic formula with an endemic intercept and population offset
 end = list(f = ~ 1, offset = population(noroBE)),
 family = "NegBin1"
 )

Now we can already fit the hhh4 model,

noroEndFit <- hhh4(stsObj = noroBE, control = noroModel)

… and directly look at the results.

summary(noroEndFit)


Call: 
hhh4(stsObj = noroBE, control = noroModel)

Coefficients:
          Estimate  Std. Error
end.1     4.18163   0.02084   
overdisp  0.89123   0.03144   

Log-likelihood:   -6883.97 
AIC:              13771.94 
BIC:              13783.57 

Number of units:        12 
Number of time points:  207

2.3 Using `RTMB`

To implement the endemic-epidemic model in RTMB we follow the introduction of the package.

The main thing is to create an objective function, that is the (log-) likelihood function that we want to maximize. As most optimizers, as well as nlminb are build for minimization, we use the negative log-likelihood.

The objective function in the simple endemic-only case can be build like this:

 objectiveEnd <- function(parms, data) {

  getAll(parms, data) # makes data list and parameters visible inside the function

  ## tells RTMB that that Y is the response
  Y <- OBS(Y)
 
  ## transform parameters
  exp_overdispersion <- exp(log_overdispersion)

  ## Endemic log linear predictor
  nu <- end.1 + log(pop[2:n_time, ])
  ## MU
  mu <- exp(nu) 
  ## negative log likelihood
  obs <- Y[2:n_time, ] # omit first time point
  nll <- - sum(dnbinom2(obs, mu = mu,  size = exp_overdispersion, log = TRUE))
  return(nll)
}

Then we need the parameter list, that also serves as inital guess when fitting the model:

parms <- list(
  end.1 = 0,
  log_overdispersion = 0
)

Since we added a data argument to the objective function objectiveEnd(), we need to build a general function to combine the objective with a specific data set. Further described in the introduction.

cmb <- function(f, d) function(p) f(p, d)

The objective function objectiveEnd() is then processed by RTMB::MakeADFun(), which builds a TMB model object compatible with the TMB package. It then constructs objective functions with derivatives from the objective function given in R (translated and done in C++) and the resulting object contains the parameters as par, functions to caluclate the objective as fn and the gradient as gr.

objEnd <- MakeADFun(cmb(objectiveEnd, data), parms)

The model is fitted using the nlminb optimizer (any gradient based optimizer in R could be used),

ptm <- proc.time() # track the runtime for optimization
optEnd <- nlminb(
  objEnd$par, # parameters
  objEnd$fn, # likelihood function
  objEnd$gr # gradient
  )
runtime_opt <- proc.time() - ptm

and uncertainties are calculated with the RTMB function sdreport().

ptm <- proc.time() # track the runtime
sdrEnd <- sdreport(objEnd)
runtime_sdr <- proc.time() - ptm

runtime_End <- runtime_opt + runtime_sdr

Note

How RTMB tracks parameters

Notice that we don’t need to manually update objEnd with the results of nlminb(). As mentioned before, MakeADFun() returns a list that contains the objective objEnd$fn and gradient functions objEnd$gr, as well as an environment objEnd$env that keeps track of the current parameter values and likelihood evaluations during the optimization performed by nlminb(). Downstream, functions like sdreport() can use objEnd immediately because it contains all that is needed.

2.4 Comparison Table

To ensure our RTMB implementation matches the hhh4() function, we can look at the parameter estimates, standard errors (SE), and runtimes in this comparison table.

Parameter estimates and runtime for <code>hhh4</code> and <code>RTMB</code> implementations with standard errors (SE) in brackets and differences between both implementations.
	hhh4	RTMB	diff.
Parameters
end.1	4.18163056 (0.02083525)	4.181630558 (0.02083525)	1.91e-09 (2.79e-10)
-log(overdisp)	0.115152025 (0.03527974)	0.115152026 (0.035279742)	-1.63e-09 (-1.24e-09)
Metrics
Log-Likelihood	-6883.968138079	-6883.968138079	0
Runtime (seconds)	0.029	0.036	-0.007

The table confirms that our RTMB implementation is mathematically equivalent to hhh4(). The parameter estimates and SEs match up to the 5th decimal and the log-likelihood is nearly the same. For this simple model formulation hhh4 is 0.007 seconds faster. Since hhh4 uses exact analytical differentiation that has been manually derived and programmed for this specific type of models, this is not suprising.

Now that the RTMB implementation of the simple endemic model is verified, we implement full endemic-epidemic models in the next section.

2 Endemic-only Models

2.1 Mathematical Formulation

2.2 Using surveillance

2.3 Using RTMB

2.4 Comparison Table

2.2 Using `surveillance`

2.3 Using `RTMB`