We describe how to partially pool parameters underlying the reproduction
numbers. This is done using a special operator in the formula passed to
epirt(). If you have previously used any of the lme4,
nlmer, gamm4, glmer or rstanarm packages then this
syntax will be familiar.
A general R formula is written as y ~ model, where y is the response that
is modeled as some function of the linear predictor which is symbolically
represented by model. model is made up of a series of terms separated by
+. In epidemia, as in many other packages, parameters can be partially pooled
by using terms of the form (expr | factor), where both expr and factor are
R expressions. expr is a standard linear model (i.e. treated the same as
model), and is parsed to produce a model matrix. The syntax (expr | factor)
makes explicit that columns in this model matrix have separate effects for
different levels of the factor variable.
Of course, separate effects can also be specified using the standard interaction
operator :. This however corresponds to no pooling, in that parameters
at different levels are given separate priors. The | operator, on the other
hand, ensures that effects for different levels are given a common prior. This
common prior itself has parameters which are given hyperpriors. This allows
information to be shared between different levels of the factor. To be concrete, suppose
that the model matrix parsed from expr has \(p\) columns,
and that factor has \(L\) levels. The \(p\)-dimensional parameter vector for
the \(l\)th group can be denoted by \(\theta_l\). In epidemia, this vector
is modeled as multivariate normal with an unknown covariance matrix. Specifically,
\[\begin{equation}
\theta_{l} \sim N(0, \Sigma),
\end{equation}\]
where the covariance \(\Sigma\) is given a prior. epidemia offers the same
priors for covariance matrices as rstanarm; in particular the decov() and lkj() priors from rstanarm can be used. Note that \(\Sigma\) is not assumed diagonal, i.e.
the effects within each level may be correlated.
If independence is desired for parameters in \(\theta_l\), we can simply replace
(expr | factor) with (expr || factor). This latter term effectively
expands into \(p\) terms of the form (expr_1 | factor), \(\ldots\), (expr_p | factor),
where expr_1 produces the first column of the model matrix given by expr,
and so on. From the above discussion, the effects are independent across terms,
and essentially \(\Sigma\) is replaced by \(p\) one-dimensional covariance matrices
(i.e. variances).
The easiest way to become familiar with how the | operator works is to
see a multitude of examples. Here, we give many examples, their interpretations,
and where possible we compare the models to the no pooling and full pooling equivalents.
For a comprehensive reference on mixed model formulas, please see Bates et al. (2015).
There are many possible ways to specify intercepts. Table 1.1 demonstrates some of these, including fully pooled, partially pooled and unpooled. Effects may also be partially pooled. This is shown in Table 1.2.
| Formula R.H.S. | Interpretation |
|---|---|
1 + ...
|
Full pooling, common intercept for all regions. |
region + ...
|
Separate intercepts for each region, not pooled. |
(1 | region) + ...
|
Separate intercepts for each region which are partially pooled. |
(1 | continent) + ...
|
Separate intercepts based on a factor other than region, partially pooled.
|
| Formula R.H.S. | Interpretation |
|---|---|
1 + npi + ...
|
Full pooling. Effect of NPI the same across all regions. |
1 + npi:region + ...
|
No pooling. Separate effect in each region. |
1 + (0 + npi|region) + ...
|
Partial pooling. Separate effects in each region. |
1 + (npi|region) + ...
|
Right hand side expands to 1 + (1 + npi|region), and so both the intercept and effect are partially pooled.
|
The final example in Table 1.2 shows that it is important to remember that to parse the term (expr | factor), epim() first parses expr into a model matrix in the same way as functions like lm() and glm() parse models. In this case, the intercept term is implicit. Therefore, if this is to be avoided, we must
explicitly use either (0 + npi | region) or (-1 + npi | region).
By default, the vector of partially pooled intercepts and slopes for each region
are correlated. The || operator can be used to specify independence. For
example, consider a formula of the form
R(region, date) ~ npi + (npi || region) + ...
The right hand side expands to 1 + npi + (1 | region) + (npi | region) + .... Separate intercepts and effects for each region which are partially pooled. The intercept and NPI effect
are assumed independent within regions.
Often groupings that are nested. For example, suppose we wish to model an epidemic at quite a fine scale, say at the level of local districts. Often there will be little data for any given district, and so no pooling will give highly variable estimates of reproduction numbers. Nonetheless, pooling at a broad scale, say at the country level may hide region specific variations.
If we have another variable, say county, which denotes the county to which
each district belongs, we can in theory use a formula of the form
R(district, date) ~ (1 | county / district) + ...
The right hand side expands to (1 | county) + (1 | county:district).
There is a county level intercept, which is partially
pooled across different counties. There are also district intercepts
which are partially pooled within each county.
Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.