# Statistically-estimated tree composition for the northeastern United States
at the time of Euro-American settlement

Christopher J. Paciorek1,*, Simon J. Goring2,&, Andrew L. Thurman3,&,
Charles V. Cogbill4, John W. Williams2,3, David J. Mladenoff6,
Jody A. Peters7, Jun Zhu8, and Jason S. McLachlan7

###### Abstract

We present a gridded 8 km-resolution data product of the estimated composition
of tree taxa at the time of Euro-American settlement of the northeastern
United States and the statistical methodology used to produce the product from
trees recorded by land surveyors. Composition is defined as the proportion of
stems larger than approximately 20 cm diameter at breast height for 22 tree
taxa, generally at the genus level. The data come from settlement-era public
survey records that are transcribed and then aggregated spatially, giving
count data. The domain is divided into two regions, eastern (Maine to Ohio)
and midwestern (Indiana to Minnesota). Public Land Survey point data in the
midwestern region (ca. 0.8-km resolution) are aggregated to a regular 8 km
grid, while data in the eastern region, from Town Proprietor Surveys, are
aggregated at the township level in irregularly-shaped local administrative
units. The product is based on a Bayesian statistical model fit to the count
data that estimates composition on a regular 8 km grid across the entire
domain. The statistical model is designed to handle data from both the regular
grid and the irregularly-shaped townships and allows us to estimate
composition at locations with no data and to smooth over noise caused by
limited counts in locations with data. Critically, the model also allows us to
quantify uncertainty in our composition estimates, making the product suitable
for applications employing data assimilation. We expect this data product to
be useful for understanding the state of vegetation in the northeastern United
States prior to large-scale Euro-American settlement. In addition to specific
regional questions, the data product can also serve as a baseline against
which to investigate how forests and ecosystems change after intensive
settlement. The data product is being made available at the NIS data portal as
version 1.0.

Keywords: biogeography, species composition, old-growth forests, spatial
modeling, Bayesian statistical model, vegetation mapping

1Department of Statistics, University of California, Berkeley, California, USA

2Department of Geography, University of Wisconsin, Madison, Wisconsin, USA

3VA Office of Rural Health, Veterans Rural Health Resource Center, Iowa City
VAMC, Iowa City, Iowa, USA

4Harvard Forest, Harvard University, Petersham, Massachusetts, USA

5Center for Climatic Research, University of Wisconsin, Madison, Wisconsin,
USA

6Department of Forest and Wildlife Ecology, University of Wisconsin, Madison,
Wisconsin, USA

7Department of Biological Sciences, University of Notre Dame, Notre Dame,
Indiana, USA

8Department of Statistics, University of Wisconsin, Madison, Wisconsin, USA

*Corresponding author; E-mail<EMAIL_ADDRESS>
&These authors contributed equally to this work.

##  1 Introduction

Historical datasets provide critical context to understand forest ecology.
They allow researchers to define ‘baseline’ conditions for conservation
management, to understand ecosystem processes at decadal and centennial
scales, to track forest responses to shifting climates, and, particularly in
regions with widespread land use change, to understand the extent to which
forests after conversion and regeneration differ from the original forest
cover.

Euro-American settlement and subsequent land use change occurred in a time-
transient fashion across North America and were accompanied by land surveys
needed to demarcate land for land tenure and use. Various systems were used by
surveyors to locate legal boundary markers, usually by recording and marking
trees adjacent to survey markers. These data provide vegetation information
that can be mapped and used quantitatively to represent the period of
settlement. Early surveys (from 1620 until 1825) in the northeastern United
States provide spatially-aggregated data at the township level (Cogbill et
al., 2002; Thompson et al., 2013), with typical township size on the order of
200 $\mbox{km}^{2}$ and no information about the locations of individual
trees; we refer to these as the Town Proprietor Survey (TPS). Later surveys
after the establishment of the U.S. Public Land Survey System (PLS) by the
General Land Office (GLO) provide point-level data along a regular grid, with
one-half mile (800 m) spacing, for Ohio and westward during the period 1785 to
1907 (Bourdo, 1956; Pattison, 1957; Schulte and Mladenoff, 2001; Goring et
al., 2016). At each point 2-4 trees were identified, and the common name,
diameter at breast height, and distance and bearing from the point were
recorded. Survey instructions during the PLS varied through time and by point
type. Accounting for this variation requires data screening to maximize
consistency among points and the application of spatially-varying correction
factors Goring et al. (2016) to accurately assess tree stem density, basal
area and biomass from the early settlement records, but the impact on
composition estimates is limited (Liu et al., 2011). Surveyors sometimes used
ambiguous common names, which requires matching to scientific names and
standardization (Mladenoff et al., 2002; Goring et al., 2016).

Logging, agriculture, and land abandonment have left an indelible mark on
forests in the northeastern United States (Foster et al., 1998; Rhemtulla et
al., 2009b; Thompson et al., 2013; Goring et al., 2016). However most studies
have assessed these effects in individual states or smaller domains (Friedman
and Reich, 2005; Rhemtulla et al., 2009a) and with various spatial
resolutions, from townships (36 square miles) to forest zones of hundreds or
thousands of square miles. Goring et al. (2016) provide a new dataset of
forest composition, biomass, and stem density based on PLS data for the upper
Midwest that is resolved to an 8 km by 8 km grid cell scale, providing broad
spatial coverage at a spatial scale that can be compared to modern forests
using Forest Inventory and Analysis products (Gray et al., 2012). Combined
with additional, coarsely-sampled PLS data from Illinois and Indiana, newly-
digitized data from southern Michigan, and with the TPS data, this gives us
raw data for much of the northeastern United States. However, there are
several limitations of using the raw data that can be alleviated by the use of
a statistical model to develop a statistically-estimated data product. First,
the PLS and TPS data only provide estimates of within-cell variance that do
not account for information from nearby locations. Second, there are data
gaps: the available digitized data from Illinois and Indiana represent a small
fraction of those states, and missing townships are common in the TPS data.
Third, the TPS and PLS data have fundamentally different sampling design and
spatial resolution. Our statistical model allows us to provide a spatially-
complete data product of settlement-era tree composition for a common 8 km
grid with uncertainty across the northeastern U.S.

In Section 2 we describe the data sources, while Section 3 describes our
statistical models. In Section 4 we quantitatively compare competing
statistical specifications, and in Section 5 we describe the final data
product. In Section 6 we discuss the uncertainties estimated by and the
limitations of the statistical model, and we list related data products under
development.

##  2 Data

The raw data were obtained from land division survey records collated and
digitized from across the northeastern U.S. by a number of researchers (Fig.
1). For the states of Minnesota, Wisconsin, Illinois, Indiana, and Michigan
(the midwestern subdomain), digitized data are available at PLS survey point
locations and have been aggregated to a regular 8 km grid in the Albers
projection. (Note that for Indiana and Illinois, at the moment trees are
associated with township centroids and then assigned to 8 km grid cells based
on the centroid, but in the near future we will have point locations available
for each tree.) For the states of Ohio, Pennsylvania, New Jersey, New York and
the six New England states (the eastern subdomain), data are aggregated at the
township level. We make predictions for all of the states listed above; these
constitute our core domain. There are also data from a single township in
Quebec and a single township in northern Delaware; these data help inform
predictions in nearby locations within our core domain, but predictions are
not made for Quebec or Delaware. Digitization of PLS data in Minnesota,
Wisconsin and Michigan is essentially complete, with PLS data for nearly all 8
km grid cells, but data in Illinois and Indiana represent a sample of the full
set of grid cells, with survey record transcription ongoing. Data for the
eastern states are available for a subset of the full set of townships
covering the domain; the TPS data for some townships were lost, incomplete, or
have not been located (Cogbill et al., 2002).

Figure 1: Spatial domain of the northeastern United States, with locations
with data shown in gray. Locations are grid cells in midwestern portion and
townships in eastern portion. In addition to locations without data being
indicated in white, grid cells completely covered in water are white (e.g., a
few locations in the northwestern portion of the domain in the states of
Minnesota and Wisconsin).

Note that surveys occurred over a period of more than 200 years as European
colonists (before U.S. independence) and the United States settled what is now
the northeastern and midwestern United States. Our estimates are for the
period of settlement represented by the survey data and therefore are time-
transgressive; they do not represent any single point in time across the
domain, but rather the state of the landscape at the time just prior to
widespread Euro-American settlement and land use (Whitney, 1996; Cogbill et
al., 2002). These forest composition datasets do include the effects of Native
American land use and early Euro-American settlement activities (e.g., Black
et al., 2006), but it is likely that the imprint of this earlier land use is
highly concentrated rather than spatially extensive (Munoz et al., 2014).

Extensive details on the upper Midwest (Minnesota, Wisconsin, Michigan) data
and processing steps are available Goring et al. (2016); key elements include
the use of only corner points, the use of only the two closest trees at each
corner point, spatially-varying correction factors for sampling effort, and a
standardized taxonomy table. The lower Midwest (Illinois, Indiana) data were
purchased from the Indiana State Archives (Indiana) and Hubtack Document
Resources (hubtack.com; Illinois) and processed using similar steps as for the
upper Midwest data. Digitization of the Illinois and Indiana data is still
underway, so many grid cells contained no data at the time the statistical
model was fit. Note that the number of trees per grid cell varies depending on
the number of survey points in a cell, with an average of 124 trees per cell.
The gridded data at the 8 km resolution for the midwest subdomain are
available through the NIS data portal (Goring and University of Wisconsin,
2016). The TPS data were compiled by C.V. Cogbill from a myriad of archival
sources representing land division surveys conducted in connection with local
settlement and are available through the NIS data portal (Cogbill, 2016a, b).

The aggregation into taxonomic groups is primarily at the genus level but is
at the species level in some cases of monospecific genera. We model the
following 22 taxa plus an “other hardwood” category: Atlantic white cedar
(_Chamaecyparis thyoides_), Ash (_Fraxinus spp._), Basswood (_Tilia
americana_), Beech (_Fagus grandifolia)_ , Birch (_Betula spp._), Black
gum/sweet gum (_Nyssa sylvatica_ and _Liquidambar styraciflua_), Cedar/juniper
(_Juniperus virginiana_ and _Thuja occidentalis_), Cherry (_Prunus spp._),
Chestnut (_Castanea dentata_), Dogwood (_Cornus spp._), Elm (_Ulmus spp._),
Fir (_Abies balsamea_), Hemlock (_Tsuga canadensis_), Hickory (_Carya spp._),
Ironwood (_Carpinus caroliniana_ and _Ostrya virginiana_), Maple (_Acer
spp._), Oak (_Quercus spp._), Pine (_Pinus spp._), Poplar/tulip poplar
(_Populus spp._ and _Liriodendron tulipifera_), Spruce (_Picea spp._),
Tamarack (_Larix laricina_), and Walnut (_Juglans spp._). Note that in several
cases (black gum/sweet gum, ironwood, poplar/tulip poplar, cedar/juniper),
because of ambiguity in the common tree names used by surveyors, a group
represents trees from different genera or even families. For the midwestern
subdomain we do not fit statistical models for Atlantic white cedar and
chestnut as these species have 0 and 7 trees present, respectively. The taxa
grouped into the other hardwood category are those for which fewer than
roughly 2000 trees were present in the dataset; however, we include Atlantic
white cedar explicitly despite it only having 336 trees in the dataset because
of specific ecological interest in Atlantic white cedar wetlands.

Diameters are only recorded in the PLS data. Although surveyors avoided using
small trees, there was no consistent lower diameter limit. The PLS data
generally represent trees greater than 8 inches (~20 cm) diameter at breast
height (dbh), but with some trees as small as 1 inch dbh (smaller trees were
much more common in far northern Minnesota). TPS data have no information
about dbh, but the trees were large enough to blaze and are presumed to be
relatively large trees useful for marking property boundaries.

There are approximately 860,000 trees in the midwestern subdomain and 420,000
trees in the eastern subdomain. In the midwestern subdomain, oak is the most
common taxon and pine the second most common, while in the eastern subdomain
oak is the most common and beech the second most common.

Our domain is a rectangle covering all of the states using a metric Albers
(Great Lakes and St. Lawrence) projection (PROJ4: EPSG:3175), with the
rectangle split into 8 km cells, arranged in a 296 by 180 grid of cells, with
the centroid of the cell in the southwest corner located at (-71000 m, 58000
m). For the midwestern subdomain we use the western-most 146 by 180 grid of
cells when fitting the statistical models. For the eastern subdomain we use
the eastern-most 180 by 180 grid of cells and then omit 23 rows of cells in
the north and 17 rows of cells in the south, as these grid cells are outside
of the states containing data.

##  3 Statistical model

We fit a Bayesian statistical model to the data, with two primary goals:

  1. 1.

To estimate composition on a regular grid across the entire domain, filling
gaps where no data are available, and

  2. 2.

To quantify uncertainty in composition at all locations. Even in grid cells
and townships with data, we wish to quantify uncertainty because the empirical
proportions represent estimates of the true proportions that could be
calculated using the full population of all the trees in a grid cell or
township.

At a high level, the Bayesian statistical model estimates composition across
the domain, even in locations with sparse or no data, by combining the raw
composition data with the assumption that composition varies in a smooth
spatial fashion across the domain. The information in the data is quantified
by the data model, also known as the likelihood. The assumption of smoothness
is built into the model by representing the true unknown spatially-varying
composition using a statistical spatial process representation that induces
smoothing of estimates across nearby locations. This spatial process
representation is a form of prior distribution and is a function of model
parameters called hyperparameters that determine the correlation structure of
the process and are also estimated based on the data.

The result of fitting the Bayesian model via Markov chain Monte Carlo (MCMC)
is a set of representative samples from the posterior distribution for the
composition in the 23 taxonomic groupings at each of the grid cells. These
samples are the data product (described further in the Section 5) and can be
used in subsequent analyses. The mean and standard deviation of the samples
for each pair of cell and taxon represent our best estimate (i.e., prediction)
of composition and a Bayesian “standard error” quantifying the uncertainty in
the estimate.

In the remainder of this section we provide the technical specification of the
model and of the computations involved in fitting the model.

###  3.1 Data model

We start by describing the basic model for those states for which we have raw
data on the 8 km grid, and in Section 3.5 we describe the extension of the
model to accommodate data aggregated at the township level.

The statistical model treats the observations as coming from a multinomial
distribution with a (latent) vector of proportions for each grid cell,

$y_{i}\sim\mbox{Multi}(n_{i},\theta(s_{i})),$

where $y_{i}$ is the vector of counts for the $P$ taxa at the $i$th cell,
$n_{i}$ is the number of trees counted in the cell, and $\theta(s_{i})$ is the
vector of unknown proportions for those taxa at that cell. Note that we use a
standard multinomial distribution without overdispersion, because the set of
trees in the dataset is roughly uniformly sampled across the cells or
townships (Goring et al., 2016).

The proportions, $\theta_{p}(s_{i}),\,p=1,\ldots,P$, are modeled spatially by
a set of $P$ Gaussian spatial processes, one per taxon,
$\alpha_{p}(s_{i}),\,p=1,\ldots,P$. This collection of processes defines a
multivariate spatial process for composition. The $\alpha_{p}(s)$ processes
are defined on the 8 km grid,
$\alpha_{p}=\\{\alpha_{p}(s_{1}),\ldots,\alpha_{p}(s_{m})\\}$ for the $m$ grid
cells. In Section 3.4 we introduce a multinomial probit model that relates the
$\alpha_{p}(s)$ processes to the proportion processes, $\theta_{p}(s)$, via
the introduction of latent variables, with an implicit sum-to-one constraint,
$\sum_{p=1}^{P}\theta_{p}(s)=1$. A multinomial probit model is similar to
logistic regression, used for modeling a binary outcome based on an underlying
probability the outcome will occur, but generalizes to modeling a categorical
variable based on probabilities for each category.

The critical component of the statistical model is the representation of
$\alpha_{p}(s)$ as a spatial process. A spatial process is a statistical
representation that models spatially-correlated values. It provides a prior
structure that serves to smooth across noise in the observations and allows
for prediction at locations based on information from nearby locations,
including interpolation to locations with no data. Apart from the sum-to-one
constraint, the taxa are considered to be independent in the prior. We did not
want to impose any structure that ties the different taxa together, as any
correlation will likely vary across space.

In the next section, we consider two spatial models to define the structure of
the $\alpha_{p}(s)$ processes, a standard conditional autoregressive model
(Banerjee et al., 2004) and a Gaussian Markov random field (MRF) approximation
to a Gaussian process with Matérn covariance (Lindgren et al., 2011). These
models are specific statistical formulations of spatial processes that
represent spatial correlation by defining neighborhoods around each location
that are used to help inform predictions at the location.

###  3.2 Spatial process models

MRF models represent the neighborhood information by working directly with the
precision matrix (the inverse of the covariance matrix) of the values of the
spatial process, so calculation of the prior density of $\alpha_{p}$ is
computationally simple (Rue and Held, 2005). However, in situations in which
the likelihood is not normal, such as our multinomial likelihood, it can be
difficult to set up effective MCMC algorithms that are able to move in the
high-dimensional space of $\alpha_{p}$. The latent variable representation of
Section 3.4 helps to alleviate this problem. Next we describe two alternative
spatial models that we considered; in Section 4, we evaluate the models on
held-out data to choose between the two.

#### Standard conditional autoregressive models

Our first model is a standard conditional autoregressive (CAR) model;
technical details can be found in (Banerjee et al., 2004). We use a standard
form of this model that treats the four cardinal neighbors of each grid cell
as the neighbors of the grid cell. The corresponding precision matrix has
diagonal elements, $Q_{ii}$, equal to the number of neighbors for the $i$th
area (i.e., four except for cells on the boundary of the domain), while
$Q_{ik}=-1$ (the negative of a weight of one) when areas $i$ and $k$ are
neighbors and $Q_{ik}=0$ when they are not. This gives the following model for
the values of $\alpha_{p}(s_{i})$ collected as a vector across all of the grid
cells, $i=1,\ldots,m$:

$\alpha_{p}\sim\mbox{N}(0,\sigma_{p}^{2}Q^{-}).$

The use of the generalized inverse notation indicates that $Q$ is not full-
rank, but is of rank $m-1$; this gives an improper prior on an implicit
overall mean for the process values. Note that we specify an explicit mean of
zero because a non-zero mean would not be identifiable in light of the
implicit mean. This specification is called an intrinsic conditional
autoregression (ICAR) and we can write $Q=D-C$ where $C$ is the $m\times m$
adjacency matrix defining the neighborhood relation of the locations; that is
$C{}_{ik}=1$ if locations $i$ and $k$ are neighbors and zero otherwise. The
matrix $D$ is an $m\times m$ diagonal matrix containing the row sums of matrix
$C$ as the diagonal entries, $D{}_{ii}={\displaystyle\sum_{k=1}^{m}C{}_{ik}}.$

We refer to this as the _CAR model_.

#### Gaussian process approximation

Gaussian processes (GP) are also standard models for spatial processes
(Banerjee et al., 2004). GP models are computationally challenging for large
datasets because of computational manipulations involving large covariance
matrices. Given this, Lindgren et al. (2011) proposed a new framework for
using Gaussian MRFs (GMRFs) as approximations to GPs, based on the use of
stochastic partial differential equations (SPDEs).

Gaussian processes are generally constructed using one of a number of
correlation functions that define how the strength of correlation between the
values of the process at two locations decays as a function of the distance
between the locations. We consider Gaussian processes in the commonly-used
Matérn class, using the following parameterization of the Matérn correlation
function,

$R(d)=\frac{1}{\Gamma(\nu)2^{\nu-1}}\left(\frac{2\sqrt{\nu}d}{\rho}\right)^{\nu}\mathcal{K}_{\nu}\left(\frac{2\sqrt{\nu}d}{\rho}\right),$

where $d$ is Euclidean distance, $\rho$ is the spatial range parameter, and
$\mathcal{K}_{\nu}(\cdot)$ is the modified Bessel function of the second kind,
whose order is the smoothness (differentiability) parameter, $\nu>0$.
$\nu=0.5$ gives the exponential covariance. For any pair of locations, $R(d)$
defines the correlation of the process, (i.e., $\alpha_{p}(s)$ in our
context), as a function of the distance between the locations. Considering all
pairs of locations, this defines a correlation matrix for all locations of
interest.

The approach of Lindgren et al. (2011) allows us to consider MRF
approximations to the Matérn -based GP for $\nu=1$ and $\nu=2$. Our second
spatial model is this Lindgren approximation for Matérn -based GPs with
$\nu=1$. To implement this Lindgren model, one modifies the $Q$ matrix defined
previously as follows based on the technical specification of the precision
matrix provided in Lindgren et al. (2011). Let $a=4+\frac{1}{\rho^{2}}$. The
diagonal elements of $Q$ are $4+a^{2}$. The entries corresponding to cardinal
neighbors are $-2a$. Those for diagonal neighbors are $2$, and those for 2nd-
order cardinal neighbors are $1$. This extends the neighborhood structure
relative to the CAR model and parameterizes it as a function of $\rho$.

The primary difference between the CAR and Lindgren models is that the
Lindgren model provides an additional degree of freedom by estimating $\rho$.
In particular $\rho$ allows us to estimate the locality of the spatial
smoothing. As $\rho$ decreases, the model uses increasingly localized data to
estimate the compositional proportions at a given location, effectively
averaging the empirical proportions over smaller neighborhoods. In general,
the Lindgren et al. (2011) model will generally provide for a smoother
estimate than the CAR model (Paciorek, 2013).

To ensure that the $\sigma^{2}$ parameter is mathematically equivalent between
the two models, we reparameterize, producing our second model:

$\alpha_{p}\sim\mbox{N}\left(\mu_{p},\sigma_{p}^{2}\cdot\frac{4\pi}{\rho_{p}^{2}}Q(\rho_{p})^{-1}\right)$

We refer to this model as the _SPDE model_.

###  3.3 Prior Distributions

The ICAR specification contains a set of hyperparameters
$\\{\sigma_{p}^{2}\\},$ $p=1,\ldots P$. Following Gelman (2006) we use a
uniform distribution on each $\sigma_{p}$ parameter, with upper bound of 1000.
For the SPDE model we also have hyperparameters $\\{\mu_{p}\\}$, which we give
flat, non-informative priors (truncated at $\pm 10$), and $\\{\rho_{p}\\}$,
which we give uniform priors on the interval $(0.1,\exp(5))$. These various
hyperparameters are unknown parameters that control the spatial structure of
the two spatial models and are estimated from both the data and the prior
distributions just specified based on the Bayesian approach.

###  3.4 Latent Variable Model

It is well-known that devising an effective MCMC algorithm for models with
latent Gaussian process(es) and a non-Gaussian likelihood is difficult (Rue
and Held, 2005; Christensen et al., 2006; Tan and Nott, 2013). To develop an
algorithm, we make use of a latent variable representation for the multinomial
probit model (McCulloch and Rossi, 1994). The representation introduces latent
variables that allow one to develop an MCMC sampling strategy that takes
advantage of closed-form full conditional distributions (so-called Gibbs
sampling steps) for $\alpha_{p}$.

Suppose that compositional counts are available at a number of locations. At
location $i$, a sample of size $n_{i}$ observations is collected, and each
observation (i.e., each tree) can be classified into $P$ distinct categories.
For a given tree $j$ at location $i$, let $Y_{ij}$ denote the response
variable indicating the category. Let $Y_{ij}$ be associated with $P$ latent
variables $W_{ij1},...,W_{ijP}$ such that $Y_{ij}$ = $p$ if and only if
$W_{ijp}={\displaystyle\max_{p^{\prime}}\big{\\{}W_{ijp^{\prime}}\big{\\}}}$;
in other words, the maximum of the set of latent variables
$\\{W_{ijp}\\}{{}_{p=1}^{P}}$ determines the category of observation $j$ at
location $i$. The final piece of the latent variable representation is the
relationship between the $W$ variables and the $\alpha_{p}(s)$ processes. We
have that

$W_{ijp}\sim\mbox{N}(\alpha_{p}(s_{i}),1)$

independently for all of the $W_{ijp}$ values. Consider the following example
with two locations that are neighbors and $P=2$ categories. Each tree $j$ at
location $i$ is associated with two variables $W_{ij1}$ and $W_{ij2}$,
governed by the latent variables $\alpha_{1}(s_{i})$ and $\alpha_{2}(s_{i})$,
respectively. Suppose that $\alpha_{1}(s_{i})>\alpha_{2}(s_{i})$ for a given
location $i$. Then this model implies that any tree $j$ is more likely to be
labeled 1 than 2 at location $i$. The difference between $\alpha_{1}(s_{i})$
and $\alpha_{2}(s_{i})$ explains the difference in probability of categories 1
and 2 at location $i$, and the similarity between $\alpha_{p}(s_{1})$ and
$\alpha_{p}(s_{2})$ explains the correlation between the probabilities at
locations 1 and 2 for category $p$.

###  3.5 Model for township data

We developed an extension of the model described in previous sections to
account for data at a different aggregation than our core 8 km grid. This
extension introduces a new set of latent variables, one per tree, that
indicate the grid cells in which the trees are located and that can be sampled
within the MCMC as additional unknown parameters. Specifically, $c_{tj}$ is
the latent “membership” variable for tree $j$ in township $t$, $t=1,\ldots,T$.
The prior for $c_{tj}$ is a discrete distribution that puts mass, $\psi_{ti}$,
$i=1,\ldots,m$, proportional to the areal overlap between the township in
which the tree is located and the $m$ grid cells, giving

$c_{tj}\sim\mbox{Multinom}(1,\\{\psi_{t1},\ldots,\psi_{tm}\\}),$

independently across all trees. Because the townships overlap a limited number
of grid cells, most of the $\psi_{t1},\ldots,\psi_{tm}$ values are zero.

Using the latent variable representation, we have that
$W_{tjp}\sim\mbox{N}(\alpha_{p}(s_{c_{tj}}),1)$ for tree $j$ in township $t$.
In updating the other parameters in the model during the MCMC (specifically
the $\alpha$ values), we condition on the current values, $\\{c_{tj}\\}$,
which provides a “soft” (i.e., probabilistic) assignment of trees to grid
cells that respects both the known township in which the trees occurred and
the uncertainty in which grid cells the trees occurred.

Note that this prior represents the location of each tree in a township as
being independent of the other trees; this is somewhat unrealistic because it
does not represent our knowledge that the trees in a township would be
distributed somewhat regularly across the area of the township because the
witness trees were used to indicate property boundaries.

###  3.6 Computation

The McCulloch and Rossi (1994) representation is convenient for MCMC sampling,
particularly in this high-dimensional spatial context, as it allows us to draw
from the posterior conditional distributions of the $W_{ijp}$ variables (these
distributions are truncated normal) in closed form and to draw the entire
vector of latent process values for each taxon, $\alpha_{p}$, as a single
sample that respects the spatial dependence structure for each taxon.

While the latent variable representation provides great advantages in the MCMC
sampling for each $\alpha_{p}$ compared to joint Metropolis updates or
updating each location individually, there is still strong dependence between
the hyperparameters, $\\{\sigma_{p}^{2},\mu_{p},\rho_{p}\\}$ and the latent
process values (as well as between the latent process values and the latent
$W_{ijp}$ variables). To address the first, we developed a “cross-level” joint
updating strategy for the CAR model in which we propose
$\phi_{p}=\sigma_{p},p=1,\ldots,P,$ (and for the SPDE model,
$\phi_{p}\in\\{\mu_{p},(\sigma_{p},\rho_{p})\\}$) via a Metropolis-style
random walk and then given the proposed value, $\phi_{p}^{*}$, propose
$\alpha_{p}$ from its full conditional distribution given $\phi_{p}^{*}$ and
the latent $W_{p}$ variables, where $W_{p}$ is the vector of all $W_{ijp}$
values for taxon $p$: $W_{p}=\\{W_{ijp}\\},i=1,\ldots,m;j=1,\ldots,n_{i}$.
This is equivalent to sampling from the marginalized (with respect to
$\alpha_{p}$) distribution of $\phi_{p}$ conditional on $W_{p}$. For these
various joint samples of hyperparameters and $\alpha_{p}$, we use adaptive
Metropolis sampling (Shaby and Wells, 2011).

The full description of the MCMC sampling steps is provided in the Appendix.
In addition, in the latent variable representation, $\theta_{p}(s)$ never
appears explicitly and cannot be calculated in close form. Instead we use
Monte Carlo integration over $W_{ijp},\,p=1,\ldots,P$ to estimate
$\theta_{p}(s_{i})$, also described in the Appendix.

The model is implemented in R (R Core Team, 2014) with core computational
calculations coded in C++ using the _Rcpp_ package (Eddelbuettel and Francois,
2011). We also make extensive use of sparse matrix representations and
algorithms, using the _spam_ package in R (Furrer and Sain, 2010). All code is
available on Github, including pre- and post-processing code, at
https://github.com/PalEON-Project/composition.

##  4 Model comparison

###  4.1 Design

We compared the CAR and SPDE models by holding out data from the fitting
process and assessing the fit of the model on the held-out data. We used two
experimentswith held-out data:

  1. 1.

The first experiment used a subregion containing most of Minnesota and a small
amount of western Wisconsin, defined to be the cells whose x-coordinate was
less than 300,000 m (this defines a north-south line that approximately goes
through Duluth, Minnesota) and hereafter referred to as the “Minnesota
subregion”. We chose this subregion for evaluation because of its high data
density, allowing us to experiment with the effects of increasing data
sparsity on model performance. We held out all data from 95% of the cells in
this Minnesota subregion, with cells selected at random. This was meant to
assess the ability of the model to interpolate from a sparse set of
cells/townships and mimics the limited data in Illinois and Indiana.

  2. 2.

We held out 5% of the trees from all of the trees in the dataset for the
midwestern subdomain (leaving aside the held-out Minnesota subregion cells).
This was meant to assess the ability of the model to estimate the composition
in cells in which data were available.

Finally, in a separate sensitivity analysis we instead left out 80% of the
cells in Minnesota subregion at random. This variation on the first experiment
above was meant to indicate whether our model comparison conclusions would be
robust as the digitization process for Illinois and Indiana progresses and
provides us with increasingly dense data.

There has been extensive work in the statistical literature on good metrics to
use to compare the predictive ability of models; these metrics are referred to
as scoring rules. A general conclusion from this work is that predictive
distributions should maximize sharpness subject to calibration. That is, the
predictive distribution should be as narrow as possible while being calibrated
such that the observations are consistent with the distribution (Gneiting et
al., 2007). When thinking in terms of prediction intervals as summaries of the
predictive distribution, we seek intervals that are as narrow as possible
while still covering the truth the expected proportion (e.g., 95% for a 95%
prediction interval) of the time.

Following the suggestions in Gneiting et al. (2007), we considered the
following metrics: Brier score, log predictive density, mean square prediction
error, mean absolute error, and coverage and length of prediction intervals.
Further details on each are given below. For experiment 1, we define
$Y_{i}=\\{Y_{i1},\ldots,Y_{iP}\\}$ as the count of all trees in held-out cell
$i$ and for experiment 2, $Y_{i}$ is the count of held-out individual trees in
the cell, while $y_{ijp}$ is an indicator variable taking value either 0 or 1
depending on whether the $j$th held-out tree in the $i$th cell is of taxon
$p$. $\hat{\theta}_{ip}=Y_{ip}/n_{i}$ is the empirical proportion in category
$p$ for the $n_{i}$ held-out trees in cell $i$. We calculated each of the
metrics in two ways. First, we used the posterior mean composition estimates
(as an evaluation of our core predictions), with $\tilde{\theta}_{p}(s)$ being
the posterior mean. Second, we averaged the metric over the posterior samples
(as an evaluation of our full data product, including uncertainty), taking
$\tilde{\theta}_{p}(s)$ to be an individual MCMC sample and then averaging the
metric over all the posterior samples.

  1. 1.

Brier score: Gneiting et al. (2007) suggest this metric, which has been in use
for decades. For multi-category as opposed to binary outcomes, this takes the
form

$\frac{1}{n}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\sum_{p=1}^{P}(y_{ijp}-\tilde{\theta}_{p}(s_{i}))^{2}$

where $n=\sum_{i=1}^{m}n_{i}$ is the total number of held-out trees for a
given experiment and $j$ indexes across held-out trees in cell $i$.

  2. 2.

Log predictive density: This metric takes the log of the probability density
of held-out observations under the fitted model,
$Y_{i}\sim\mbox{Multinom}(n_{i},\\{\tilde{\theta}_{1}(s_{i}),\ldots,\tilde{\theta}_{P}(s_{i})\\})$,
summing on the log scale across all of the held-out data.

While in principle, this metric should be optimal (Krnjajić and Draper, 2014),
it is very sensitive to small predictions near zero (Gneiting et al., 2007).
Even worse, our Monte Carlo estimation of $\theta$ used 10000 samples, so in
some cases $\tilde{\theta}_{p}(s)=0$. When a tree is present in a cell but its
corresponding proportion is 0, this gives a log density of $-\infty$,
preventing use of the metric. As an informal solution to this we set
$\tilde{\theta}_{p}(s)=\frac{1}{100000}$ in such cases, but given these issues
we treat the log predictive density as a secondary measure.

  3. 3.

(Experiment 1 only) Weighted root mean square prediction error (RMSPE),

$\sqrt{\frac{1}{Pn}\sum_{i=1}^{m}\sum_{p=1}^{P}n_{i}(\hat{\theta}_{ip}-\tilde{\theta}_{p}(s))^{2}}$

and mean absolute error (MAE)

$\frac{1}{Pn}\sum_{i=1}^{m}\sum_{p=1}^{P}n_{i}|\hat{\theta}_{ip}-\tilde{\theta}_{p}(s)|:$

These metrics calculate the error of the estimated proportions relative to the
empirical proportions based on the held-out trees, averaging over cells and
taxa. We weight by the number of held-out trees in each cell to account for
the greater variability in the empirical proportions in locations with few
held-out trees.

  4. 4.

(Experiment 1 only) Coverage and length of 95% prediction intervals for
$Y_{ip}$. We considered only cells with at least 50 trees to focus our
assessment on cases where empirical proportions were reasonably certain and
avoid being strongly influenced by predictive inference for cells where
observational variability dominates.

Note that all of the metrics except coverage and interval length can be
applied to individual posterior samples and therefore allow us to estimate the
posterior probability that one model has a lower (better) value of the metric
than the other model by simply calculating the proportion of samples for which
the model has a lower value of the metric. Also note that in addition to
allowing comparison between models the MAE and RMSPE metrics allow one to
assess absolute performance of each model in predicting composition.

In our initial exploratory fitting, we noticed that the SPDE model produced
boundary effects in the predicted composition near the edges of the convex
hull of the observations. To attempt to alleviate this, we added a buffer zone
with a width of six grid cells around our entire original domain, but note
that the boundary effects were still evident even after inclusion of the
buffer. For the model comparison, we included this buffer for both the SPDE
and CAR models.

We ran each model for 150,000 iterations. After discarding 25,000 iterations
for burn-in, we retained a posterior sample of 250 subsampled iterations – we
use a subsample instead of the full 125,000 post-burn-in iterations to reduce
post-processing computations and storage needs.

###  4.2 Results

Here we summarize the results of our analyses that inform the choice between
the CAR and SPDE models.

#### Full cell hold-out experiment

For Experiment 1 (full cells held out) for cells in the Minnesota subregion
held out of the fitting process, the CAR model outperforms the SPDE model
based on the posterior distribution over the predictive metric values (Table
1). For the posterior mean predictions, the SPDE model appears to outperform
the CAR model to a lesser degree, but we do not have any uncertainty estimates
for this comparison. Coverage and interval lengths are similar between the two
models (Table 2). From a practical perspective, based on the difference in
mean absolute error, the differences between the models are small (Table 1).

Table 1: Predictive ability based on several predictive metric criteria for the CAR and SPDE spatial models when holding out 95% of entire cells of data in Minnesota. Smaller values are better. | Posterior mean of metric | Metric of posterior mean predictions
---|---|---
| CAR model | SPDE model |  Posterior Prob. CAR < SPDE |  CAR model |  SPDE model
Brier | 0.819 | 0.844 |  0.98 |  0.738 |  0.733
Negative Log Density | 466325 | 510383 |  1.00 |  394003 |  394554
Mean Absolute Error | 0.0364 | 0.0383 |  0.98 |  0.0275 |  0.0269
Root Mean Square Error | 0.0897 | 0.0960 |  0.97 |  0.0647 |  0.0627
Table 2: Coverage and length of prediction intervals for the CAR and SPDE spatial models when holding out 95% of entire cells of data in Minnesota. Coverage values near 0.95 are optimal, while shorter intervals are better. | CAR model | SPDE model
---|---|---
Coverage | 0.977 | 0.978
Mean Interval Length | 0.129 | 0.142
Median Interval Length | 0.037 | 0.033

The results for the variation on Experiment 1 in which the proportion of cells
that are held out decreases from 95% to 80% show that the SPDE model generally
outperforms the CAR model, but again differences from a practical perspective,
based on mean absolute error, are limited (Tables 3-4).

Table 3: Predictive ability based on several predictive metric criteria for the CAR and SPDE spatial models when holding out 80% of entire cells of data in Minnesota. Smaller values are better. | Posterior mean of score | Score of posterior mean predictions
---|---|---
| CAR model | SPDE model |  Posterior Prob. CAR < SPDE | CAR model | SPDE model
Brier | 0.773 | 0.765 |  0.10 | 0.710 | 0.710
Negative Log Density | 355928 | 353987 |  0.25 | 311525 | 311902
Mean Absolute Error | 0.0309 | 0.0296 |  0.10 | 0.0226 | 0.0223
Root Mean Square Error | 0.0763 | 0.0739 |  0.02 | 0.0533 | 0.0530
Table 4: Coverage and length of prediction intervals for the CAR and SPDE spatial models when holding out 80% of entire cells of data in Minnesota. Coverage values near 0.95 are optimal, while shorter intervals are better. | CAR model | SPDE model
---|---|---
Coverage | 0.981 | 0.972
Mean Interval Length | 0.112 | 0.103
Median Interval Length | 0.028 | 0.022

#### Individual tree hold-out experiment

In Experiment 2 (individual trees held out), we have evidence (posterior
probability of 0.93) that the SPDE model is better based on the Brier score,
but the Brier score values for the two models are numerically almost the same
(Table 5).

Table 5: Predictive ability based on several predictive metric criteria for the CAR and SPDE spatial models when holding out 5% of trees. Smaller values are better. | Posterior mean of metric | Metric of posterior mean predictions
---|---|---
| CAR model | SPDE model |  Posterior Prob. CAR < SPDE | CAR model | SPDE model
Brier | 0.662 | 0.661 |  0.07 | 0.657 | 0.657
Negative Log Density | 51757 | 51626 |  0.01 | 50705 | 50736

#### Choice of spatial model

The differences between models are not consistent across the various
comparisons, so there is not a clear choice. In our final data product we use
the CAR model, for three reasons. First, the CAR model has modestly better
performance when data are sparse, as is still the case for Illinois and
Indiana. Second, the model is simpler and easier to explain, and computations
can be done more quickly. Third, predictions from the SPDE model showed
boundary effects, with some taxa showing non-negligible posterior mean values
at the edges of the domain, well away from where the taxa were present in the
empirical data. This included non-negligible values within (but near the edge
of) the convex hull of locations with data.

##  5 Data product

The final data product is a dataset that contains 250 posterior samples of the
proportions of each of the 23 tree taxa at each grid cell in the states in our
domain of the northeastern United States.

For this final data product, we ran the model using the CAR specification with
all of the data (including the data held out in the model comparison analyses)
for 150,000 iterations with the same burn-in and subsampling details as
described in Section 4. Based on graphical checks and calculation of effective
sample size values, mixing was generally reasonable, but for some of the
hyperparameters was relatively slow, particularly for less common taxa.
Despite this, mixing for the variables of substantive interest – the
proportions – was good, with effective sample sizes for the final product
generally near 250.

Maps of estimated composition for the full domain for several taxa of
substantive interest illustrate the results, contrasting the raw data
proportions, the posterior means, and posterior standard deviations as
pointwise estimates of uncertainty (Fig. 2). We also present the posterior
means for all 23 taxa (Fig. 3).

Figure 2: Empirical proportions from raw data (column 1), predictions in the
form of posterior means (column 2) and uncertainty estimates in the form of
posterior standard deviations – representing standard errors of prediction
(column 3) for select taxa.  Figure 3: Predictions (posterior means) for all
taxa over the entire domain.

The data product is publicly available at the NIS Data Portal under the CC BY
4.0 license as version 1.0 as of January 2016 (Paciorek et al., 2016). The
product is in the form of a netCDF-4 file, with dimensions x-coordinate,
y-coordinate, and MCMC iteration. There is one variable per taxon. In
addition, dynamic visualizations of the product using the Shiny tool are
available at https://www3.nd.edu/$\sim$paleolab/paleonproject/maps. The PalEON
Project (in particular the first author) will continue to maintain this
product, releasing new versions as additional data in Illinois, Indiana and
Ohio are digitized. Note that digitization of data from Illinois and Indiana
is ongoing, and digitization of additional data from Ohio is planned as well.
As a result, at some point we expect to have complete data for the midwestern
half of the domain.

##  6 Discussion

In the parts of the modeled region with spatially complete data (in
particular, Minnesota, Wisconsin, and Michigan), the statistical estimates of
forest composition closely match the patterns apparent in the raw data (Fig.
2), as expected. In these areas, the estimated tree composition from the model
has the advantage of downweighting unusual or outlier values in the empirical
proportions of individual grid cells, which are likely due to stochastic
sampling variation within that grid cell (compare the first two columns in
Fig. 2). Some stochastic variation is expected given that, even in the most
spatially complete regions, each grid cell contains an average of 124 trees
(120-135 is typical) (Goring et al., 2016) and some cells contain many fewer
trees. Hence, some smoothing of this stochastic variation is appropriate. This
smoothing is based on information on data from nearby cells, and the estimates
from the model reflect the smooth trends in forest composition across the
spatial domain. A partial cost is that these maps can smooth out sharp
ecotones or other forms of true spatial heterogeneity, particularly in areas
with sparse data (including areas with low tree density). For example, the
sharp increase in Elm along the Minnesota River (Fig. 2, first column) likely
represents a real ecotone in the settlement-era vegetation. Vegetation
gradients and ecotones were sharper in the settlement-era forests in the upper
Midwest than in contemporary forests (Goring et al., 2016), and the modeled
estimates may partially obscure this change. Users interested in using the
original unsmoothed data are directed to the data product described earlier
(Goring and University of Wisconsin, 2016). Additional investigation of other
statistical representations to better capture sharp gradients is of interest,
in particular nonstationary spatial models and use of covariates. Potential
environmental covariates include soils, firebreaks, and topography (Grimm,
1984; Shea et al., 2014). Here, however, we chose to limit our model to be
exclusively a function of spatial distance in order to avoid dependence on the
environmental drivers of pre-settlement forest composition that might lead to
circular reasoning in subsequent inferences based on this dataset. Use of
covariates could also lead to prediction that a taxa is present well beyond
its range boundary in places where data are sparse.

A key advance of this work over prior reconstructions of settlement-era
vegetation lies in the estimates of uncertainty across the spatial domain.
These estimates of uncertainty include the sampling uncertainty within grid
cells (as do the within-grid cell estimates of uncertainty available from the
raw proportions), but, because this is a spatial model, predictions and their
associated uncertainty estimates are also informed by the information content
of nearby cells. The maps of standard errors across species (Fig. 2, third
column) highlight the advantages of this approach in areas of high data
coverage (Minnesota, Wisconsin, Michigan) and in areas of sparse coverage
(e.g., Illinois, Indiana, parts of Ohio). Where there are not large gaps in
the data, the model provides low and fairly smooth estimates of uncertainty.
Uncertainty is generally higher in the eastern subdomain than in the areas of
the midwestern subdomain with high data coverage because of missing townships
and lower sampling density even in townships with data. In areas of sparse
coverage and in areas with low tree density (e.g., southwestern Minnesota),
the standard error of our estimates increases appropriately. Nevetheless,
these uncertainties surround reasonable estimates of trends in composition.
For example, the model does a good job of capturing the oak ecotone in Indiana
and Illinois, representing a shift from oak savannas and woodlands to closed
mesic forests (Fig. 2). Experiment 1 showed that both models predicted
composition at cells with no data reasonably well, mimicking the case of
sparsely sampled data and giving confidence in the broad spatial patterns
predicted in more poorly sampled regions, particularly those with regular, but
sparse sampling that mimic the experiment (Illinois and Indiana, but not
Ohio). The apparent blockiness of uncertainty estimates in a few places such
as Ohio is caused by spatial gaps and variations in sampling resolution.
Absolute uncertainty generally increases with abundance for all taxa (Fig. 2,
column 3).

The exploration of alternative approaches to spatial modeling of composition
showed similar results for the SPDE and CAR models, both in terms of
prediction accuracy and performance of prediction intervals. Small differences
among the various metrics of goodness of fit favored each model in turn, but
applied users of the models would find little pragmatic difference affecting
scientific inference. Ultimately, we slightly favor the CAR model, because it
avoids the boundary effects apparent in the SPDE model at the edges of the
domain.

The models presented here estimate only the relative abundance of tree taxa,
which does not directly tell us about tree density or other aspects of
vegetation structure. This becomes a particular limitation for interpreting
vegetation where trees become sparse at the prairie-forest transition from
northern Minnesota through southern Illinois (Transeau, 1935). Our model
(correctly) predicts that tree composition there is dominated by oak, but this
is less useful considering the sparseness of trees. This limitation can be
addressed by developing estimates of absolute abundance (e.g., biomass) rather
than compositional estimates. A gridded dataset of biomass, stem density, and
basal area is already available for Minnesota, Wisconsin, and northern
Michigan (Goring et al., 2016), based on the PLS data. An extension to
southern Michigan, Illinois, and Indiana is planned. We are currently
developing statistical estimates of biomass for Minnesota, Wisconsin, and
Michigan using a statistical model applied to the gridded biomass dataset,
with extension to Illinois and Indiana planned. We also plan to estimate stem
density and basal area using a similar approach to that used for biomass.

## Author contributions

CJP, ALT, and JZ developed the statistical model and code. CJP carried out the
model comparison and created the data product. CJP wrote the paper based on an
initial draft by ALT and with feedback and editing from SJG, CVC, JWW, DJM,
JAP, JZ and JSM. SJG, JAP, CVC, DJM, JSM, and JWW led the processing and
analysis of the PLS and TPS data and assisted with interpretation of results.

## Acknowledgments

The authors are deeply indebted to all of the researchers over the years who
have preserved, collected, and digitized survey records, in particular John
Burk, Jim Dyer, Peter Marks, Robert McIntosh, Ed Schools, Ted Sickley, Ronald
Stuckey, and the Ohio Biological Survey. We thank Madeline Ruid, Benjamin
Seliger, Morgan Ripp and Daniel Handel for processing of the southern Michigan
data. Indiana and Illinois data were made possible through the hard work of
many Notre Dame undergraduates in the McLachlan lab. This work was carried out
by the PalEON Project with support from the National Science Foundation
MacroSystems Program through grants EF-1065702, EF-1065656, DEB-1241874 and
DEB-1241868 and from the Notre Dame Environmental Change Initiative.

## References

  * Banerjee et al. (2004) Banerjee, S., B. Carlin, and A. Gelfand (2004).  Hierarchical Modeling and Analysis for Spatial Data.  Boca Raton, Florida: Chapman & Hall.
  * Black et al. (2006) Black, B. A., C. M. Ruffner, and M. D. Abrams (2006).  Native American influences on the forest composition of the Allegheny Plateau, northwest Pennsylvania.  Canadian Journal of Forest Research 36(5), 1266–1275.
  * Bourdo (1956) Bourdo, E. A. (1956).  A review of the General Land Office survey and of its use in quantitative studies of former forests.  Ecology 37, 754–768.
  * Christensen et al. (2006) Christensen, O., G. Roberts, and M. Sköld (2006).  Robust Markov chain Monte Carlo methods for spatial generalized linear mixed models.  Journal of Computational and Graphical Statistics 15, 1–17.
  * Cogbill (2016a) Cogbill, C. (2016a).  Settlement-era tree composition, eastern US: Level 1.  Technical Report http://dx.doi.org/10.6073/pasta/e40c1ad172882d5113844296611451f6, Long Term Ecological Research Network.
  * Cogbill (2016b) Cogbill, C. (2016b).  Settlement-era tree composition, Ohio: level 1.  Technical Report http://dx.doi.org/10.6073/pasta/d09bb56c6af8ef0783cd7e76f8113b34, Long Term Ecological Research Network.
  * Cogbill et al. (2002) Cogbill, C., J. Burk, and G. Motzkin (2002).  The forests of presettlement New England, USA: spatial and compositional patterns based on town proprietor surveys.  Journal of Biogeography 29, 1279–1304.
  * Eddelbuettel and Francois (2011) Eddelbuettel, D. and R. Francois (2011).  Rcpp: Seamless R and C++ integration.  Journal of Statistical Software 40, 1–18.
  * Foster et al. (1998) Foster, D. R., G. Motzkin, and B. Slater (1998).  Land-use history as long-term broad-scale disturbance: regional forest dynamics in central New England.  Ecosystems 1(1), 96–119.
  * Friedman and Reich (2005) Friedman, S. K. and P. B. Reich (2005).  Regional legacies of logging: departure from presettlement forest conditions in northern Minnesota.  Ecological Applications 15(2), 726–744.
  * Furrer and Sain (2010) Furrer, R. and S. R. Sain (2010).  spam: A sparse matrix R package with emphasis on MCMC methods for Gaussian Markov random fields.  Journal of Statistical Software 36(10), 1–25.
  * Gelman (2006) Gelman, A. (2006).  Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).  Bayesian Analysis 1(3), 515–534.
  * Gneiting et al. (2007) Gneiting, T., F. Balabdaoui, and A. E. Raftery (2007).  Probabilistic forecasts, calibration and sharpness.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69(2), 243–268.
  * Goring and University of Wisconsin (2016) Goring, S. and University of Wisconsin (2016).  Settlement-era gridded tree composition, midwestern US: Level 1.  Technical Report http://dx.doi.org/10.6073/pasta/aa0ef9828d41569a96651b056ad89fb3, Long Term Ecological Research Network.
  * Goring et al. (2016) Goring, S. J., J. W. Williams, D. J. Mladenoff, C. V. Cogbill, S. Record, C. J. Paciorek, S. T. Jackson, M. C. Dietze, J. H. Matthes, and J. S. McLachlan (2016).  Changes in forest composition, stem density, and biomass from the settlement era to present in the upper Midwestern United States.  PLOS ONE in press.
  * Gray et al. (2012) Gray, A. N., T. J. Brandeis, J. D. Shaw, W. H. McWilliams, and P. D. Miles (2012).  Forest inventory and analysis database of the United States of America (FIA).  Vegetation databases for the 21st century.–Biodiversity & Ecology 4, 255–264.
  * Grimm (1984) Grimm, E. C. (1984).  Fire and other factors controlling the Big Woods vegetation of Minnesota in the mid-nineteenth century.  Ecological Monographs 54(3), 291–311.
  * Krnjajić and Draper (2014) Krnjajić, M. and D. Draper (2014).  Bayesian model comparison: Log scores and DIC.  Statistics and Probability Letters 88, 9–14.
  * Lindgren et al. (2011) Lindgren, F., H. Rue, and J. Lindström (2011).  An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society B 73, 423–498.
  * Liu et al. (2011) Liu, F., D. J. Mladenoff, N. S. Keuler, and L. S. Moore (2011).  Broadscale variability in tree data of the historical public land survey and its consequences for ecological studies.  Ecological Monographs 81(2), 259–275.
  * McCulloch and Rossi (1994) McCulloch, R. and P. E. Rossi (1994).  An exact likelihood analysis of the multinomial probit model.  Journal of Econometrics 64(1), 207–240.
  * Mladenoff et al. (2002) Mladenoff, D. J., S. E. Dahir, E. V. Nordheim, L. A. Schulte, and G. G. Guntenspergen (2002).  Narrowing historical uncertainty: Probabilistic classification of ambiguously identified tree species in historical forest survey data.  Ecosystems 5(6), 539–553.
  * Munoz et al. (2014) Munoz, S. E., D. J. Mladenoff, S. Schroeder, and J. W. Williams (2014).  Defining the spatial patterns of historical land use associated with the indigenous societies of eastern North America.  Journal of Biogeography 41(12), 2195–2210.
  * Paciorek (2013) Paciorek, C. (2013).  Spatial models for point and areal data using Markov random fields on a fine grid.  Electronic Journal of Statistics 7, 946–972.
  * Paciorek et al. (2016) Paciorek, C., S. Goring, A. Thurman, C. Cogbill, J. Williams, D. Mladenoff, J. Peters, J. Zhu, and J. McLachlan (2016).  Settlement-era gridded tree composition, northeastern US: Level 2.  Technical Report http://dx.doi.org/10.6073/pasta/8544e091b64db26fdbbbafd0699fa4f9, Long Term Ecological Research Network.
  * Pattison (1957) Pattison, W. D. (1957).  Beginnings of the American Rectangular Land Survey System, 1784-1800.  Chicago: University of Chicago.
  * R Core Team (2014) R Core Team (2014).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Rhemtulla et al. (2009a) Rhemtulla, J. M., D. J. Mladenoff, and M. K. Clayton (2009a).  Historical forest baselines reveal potential for continued carbon sequestration.  Proceedings of the National Academy of Sciences 106(15), 6082–6087.
  * Rhemtulla et al. (2009b) Rhemtulla, J. M., D. J. Mladenoff, and M. K. Clayton (2009b).  Legacies of historical land use on regional forest composition and structure in Wisconsin, USA (mid-1800s-1930s-2000s).  Ecological Applications 19(4), 1061–1078.
  * Rue and Held (2005) Rue, H. and L. Held (2005).  Gaussian Markov Random Fields: Theory and Applications.  Boca Raton: Chapman & Hall.
  * Schulte and Mladenoff (2001) Schulte, L. A. and D. J. Mladenoff (2001).  The original US public land survey records: their use and limitations in reconstructing presettlement vegetation.  Journal of Forestry 99(10), 5–10.
  * Shaby and Wells (2011) Shaby, B. and M. Wells (2011).  Exploring an adaptive Metropolis algorithm.  Technical Report 2011-14, Department of Statistics, Duke University.
  * Shea et al. (2014) Shea, M. E., L. A. Schulte, and B. J. Palik (2014).  Reconstructing vegetation past: pre-Euro-American vegetation for the midwest driftless area, USA.  Ecological Restoration 32(4), 417–433.
  * Tan and Nott (2013) Tan, L. S. and D. J. Nott (2013).  Variational inference for generalized linear mixed models using partially noncentered parametrizations.  Statistical Science 28(2), 168–188.
  * Thompson et al. (2013) Thompson, J. R., D. N. Carpenter, C. V. Cogbill, and D. R. Foster (2013).  Four centuries of change in northeastern United States forests.  PLOS ONE 8(9), e72540.
  * Transeau (1935) Transeau, E. N. (1935).  The prairie peninsula.  Ecology 16(3), 423–437.
  * Whitney (1996) Whitney, G. G. (1996).  From Coastal Wilderness to Fruited Plain: a History of Environmental Change in Temperate North America from 1500 to the Present.  Cambridge University Press.

##  7 Appendix

###  7.1 MCMC details

Define $\bar{w}_{ip}=\frac{1}{n_{i}}{\displaystyle\sum_{j=1}^{n_{i}}W_{ijp}}$
as the average of the $W$ values for the $p$th taxon in the $i$th grid cell
and $\bar{w}_{p}=\\{\bar{w}_{ip}\\},$ $i=1,\ldots,m$. Let $A$ be a diagonal
matrix where $A_{ii}$ is the number of trees in the $i$th grid cell. When
there are no trees in a grid cell, $\bar{w}_{ip}=0$ and $A_{ii}=0$. For the
township data, at each iteration, based on the current values of the grid cell
membership variables, $\\{c_{tj}\\}$, trees are aggregated into grid cells and
the calculations above can then be carried out.

The conditional distribution for $W_{ijp}$ given the other unknowns in the
model and the data is as follows. Let $\mbox{TN}(a,b,\mu,\tau^{2})$ denote the
truncated normal distribution with mean parameter $\mu$ and variance parameter
$\tau^{2}$, truncated below by $a$ and above by $b$.

$\displaystyle W_{ijp}$
$\displaystyle\sim\begin{cases}\mbox{TN}\big{(}{\displaystyle\max_{p^{*}\neq
y_{ij}}w_{ijp^{*}},\infty,\alpha_{y_{ij}}(s_{i}),1\big{)},}&\mbox{if }p\mbox{
= $y_{ij}$}\\\
\mbox{TN}\big{(}-\infty,w_{ijy_{ij}},\alpha_{p}(s_{i}),1\big{)},&\mbox{if
}p\mbox{ $\neq y_{ij}$}\end{cases}$

In essence, the truncation value is determined by the taxon of the $j$th tree.
For a given $p$, the $W$ values for all trees in all cells can be sampled in
parallel.

The conditional distribution of $\alpha_{p}$ is

$\displaystyle\alpha_{p}$
$\displaystyle\sim\mbox{N}\bigg{(}\Big{(}A+Q_{p}\Big{)}^{-1}A\bar{w}_{p},\Big{(}A+Q_{p}\Big{)}^{-1}\bigg{)}.$

where $Q_{p}=(\sigma_{p}^{2})^{-1}Q$ for the CAR model and
$\left(\sigma_{p}^{2}\cdot\frac{4\pi}{\rho_{p}^{2}}\right)^{-1}Q(\rho_{p})$
for the SPDE model. For each hyperparameter, $\phi_{p}=\log\sigma_{p}$ for the
CAR model and $\phi_{p}\in\\{\mu_{p},(\log\sigma_{p},\log\rho_{p})\\}$ for the
SPDE model, we sample $\\{\phi_{p},\alpha_{p}\\}$ jointly, proposing
$\phi_{p}$ as a random walk and, conditional on the proposed value of
$\phi_{p}$, sampling $\alpha_{p}$ from the distribution just above. The joint
proposal is accepted or rejected as a standard Metropolis-Hastings proposal,
with adaptation of the proposal (co)variance (Shaby and Wells, 2011). The
proposal distribution for $\phi_{p}$ is a normal distribution (bivariate for
$\phi_{p}=(\log\sigma_{p},\log\rho_{p})$).

For the township-level data, for a given tree $j$ in township $t$, we draw the
latent tree membership variable, $c_{tj}\in\\{1,\ldots,m\\}$, from a discrete
distribution by normalizing posterior weights,
$\\{\psi_{1}L_{tj1},\ldots,\psi_{m}L_{tjm}\\}$, produced by multiplying the
prior weights by a likelihood contribution, $L_{tji}$, $i=1,\ldots,m$.
$L_{tji}$ is the density of the latent $W_{tj1},\ldots,W_{tjP}$ values for the
given tree under the condition that $c_{tj}=i$, namely the product of
independent normal densities, $W_{tjp}\sim\mbox{N}(\alpha_{p}(s_{i}),1)$, over
$p=1,\ldots,P$. Thus the posterior reweights the prior based on how consistent
the current $W_{tj}$ values for a tree are with the $\alpha$ values for the
candidate grid cells.

###  7.2 Estimating $\theta_{p}(s)$ via Monte Carlo integration

In the latent variable representation, $\theta_{p}(s)$ never appears
explicitly and cannot be calculated in closed form. Instead we use Monte Carlo
integration over $W_{ijp},\,p=1,\ldots,P$ to estimate $\theta_{p}(s_{i})$. The
quantity
$\theta_{p}(s_{i})=\mbox{Prob}(W_{ijp}={\displaystyle\max_{p^{*}}W_{ijp^{*}})}$
defines the probability of taxon $p$ at grid cell $i$. This requires one to
choose the number of Monte Carlo samples, which we set at 10000, effectively
sampling 10000 hypothetical trees and estimating the probabilities of the
different taxa in the population from the empirical proportions in this sample
of trees. For each of the saved MCMC samples, $k=1,\ldots K$, we estimate
$\theta_{p}^{(k)}(s_{i})$ numerically. Specifically, for $t=1,\ldots,10000$
samples (i.e., hypothetical trees), we independently draw

$W_{itp}^{(k)}\sim\mbox{N}(\alpha_{p}^{(k)}(s_{i}),1),\,p=1,\ldots,P$

and estimate using

$\theta_{p}^{(k)}(s_{i})\approx\frac{1}{10000}{\displaystyle\sum_{t=1}^{10000}1(W_{itp}^{(k)}={\displaystyle\max_{p^{*}}W_{itp^{*}}^{(k)})}},\,p=1,\ldots,P$

where $1(\cdot)$ is the indicator function that evaluates to 1 if the
expression is true and 0 if false. In other words, we calculate the proportion
of times that the maximum of $W_{itp},\,p=1,\ldots,P$ corresponds to taxon
$p$. Considering $\theta_{p}^{(k)}(s_{i}),\,k=1,\ldots,K$, we have a sample
from the posterior of $\theta_{p}(s_{i})$.


# Scalable Microservice Forensics and Stability Assessment Using Variational
Autoencoders

Prakhar Sharma, Phillip Porras, Steven Cheung, James Carpenter, Vinod
Yegneswaran Computer Science Laboratory, SRI International

###### Abstract

We present a deep-learning-based approach to containerized application runtime
stability analysis, and an intelligent publishing algorithm that can
dynamically adjust the depth of process-level forensics published to a backend
incident analysis repository. The approach applies variational autoencoders
(VAEs) to learn the stable runtime patterns of container images, and then
instantiates these container-specific VAEs to implement stability detection
and adaptive forensics publishing. In performance comparisons using a
50-instance container workload, a VAE-optimized service versus a conventional
eBPF-based forensic publisher demonstrates 2 orders of magnitude (OM) CPU
performance improvement, a 3 OM reduction in network transport volume, and a 4
OM reduction in Elasticsearch storage costs. We evaluate the VAE-based
stability detection technique against two attacks, CPUMiner and HTTP-flood
attack, finding that it is effective in isolating both anomalies. We believe
this technique provides a novel approach to integrating fine-grained process
monitoring and digital-forensic services into large container ecosystems that
today simply cannot be monitored by conventional techniques.

##  I Introduction

Application containerization is becoming the dominant technology to architect
and deploy microservice applications and software services. IDC Research
predicts that by 2023 more than 500 million applications and services will be
developed using cloud-native containerization services [4], which is roughly
the same number of applications that have been developed over the last 40
years combined. Containerization has become integral to the security and
management of front-line software services and sensitive data processing
applications across all major industries in our economy.

However, container-based microservice ecosystems pose a myriad of challenges
to the integration of security services, runtime management, and forensics
collection for security incident response; most significant of which is
scaling these services to manage modern container pipelines. Process
monitoring and forensic logging impose significant costs to container workload
overhead, as they consume CPU cycles, network bandwidth, and data storage, to
deliver and store information to log repositories. In fact, industrial scale
container ecosystems (thousands of application instances and beyond) are
largely out of reach of conventional process-level monitoring and forensic
services. Nevertheless, fine-grained process monitoring remains an important
service, and indeed a compliance requirement, for enabling security, fault
management, and damage assessment.

This paper examines two critical questions pertaining to the secure runtime
management of such large-scale container ecosystems. First, how can we conduct
massively scalable monitoring of containerized applications to identify when
containers become unstable or subverted? Second, when runtime instabilities
are perceived, how can we capture enough fine-grained digital forensics to
drive an incident response or fault analysis?

For container DEVOPS environments, process visibility is commonly achieved
through some form of extended Berkeley Packet Filter (eBPF)-based tracing,
such as through services like Sysdig[7] and Chisel[19]. These tools offer a
significant performance advantage over typical host auditing services for
process system-call logging. eBPF extends a host kernel with logic to trap key
system calls on entry or exit from the kernel, providing a fine-grained view
of process-internal runtime behavior to drive intrusion detection, fault
analysis, and other forms of anomaly detection.

We propose a new strategy for conducting a computationally-efficient local
analysis of container runtime behavior. This is done through an extraction of
activity vectors that captures statistics and attributes from a container’s
eBPF-based forensics produced during intervals of runtime operation. An
activity vector is input to a VAE that determines the extent to which the
vector fits the container’s trained stable baseline runtime model. VAEs are
robust and efficient in performing event-based pattern analysis [21, 22]. We
introduce a VAE-driven adaptive forensic publisher, which produces a runtime
stability assessment to determine the depth of forensics published to the
remote log repository. During an interval of runtime stability, the VAE-based
forensic publisher delivers an activity model that captures the modality of
the application as it was observed through its eBPF-derived system-call
forensics. An activity model is a highly reduced representation of the
forensic data, requiring a fraction of transport and storage cost. During
intervals not matching the trained runtime pattern, the publisher delivers the
interval’s forensic log along with the activity model. Analytics are then
performed on the forensic logs to isolate evidence of errors or faults, as
well system-call actions that may indicate malicious operations performed
during the unstable runtime interval.

We present an overview of the VAE algorithm for interval-based container
stability assessment and our adaptive forensic publisher service. The system
has been implemented and tested within Docker and Kubernetes environments.
Results of a comparative performance assessment of the VAE-based forensic
publisher versus conventional process-level forensic collection services are
presented. These results demonstrate the potential for deep-learning optimized
forensic publishing services to enable fine-grained process monitoring and
forensic-based incident response service to scale within ecosystems that today
are beyond the resource consumption limitations of conventional monitoring
techniques.

Fig 1. Design overview of VAE-based adaptive forensic publishing service for
containerized applications.

##  II VAE Design

Given any dataset $X$ drawn i.i.d from true data distribution $p_{data}$,
variational auto-encoders [8] (VAEs) aim to learn an encoder decoder pair that
enable drawing new samples $x_{new}$ such that $x_{new}\sim p_{data}$ ($\sim$
signifies a sample drawn from a probability distribution). The encoder and
decoders are usually represented by neural networks. The generative process of
VAEs involves learning a low dimensional latent space and is defined as:

$z_{new}\sim p(Z)\ \ \ x_{new}\sim p_{\theta}(X|Z=z_{new})$

where $p(Z)$ is a fixed prior distribution over latent space $Z$. VAE
operation can be broken down into two parts, first the input is mapped to the
latent space via an encoder $E$ (parametrized by $\phi$)

$E_{\phi}(x)=z\sim q_{\phi}(z|x)=q(Z|f_{\phi}(x))$

and then the latent space is mapped to the original distribution via a decoder
$D$ (parametrized by $\theta$)

$D_{\theta}(z)=x_{new}\sim p_{\theta}(x|z)=p(X|g_{\theta}(z))$

Here $f_{\phi}(x)$ and $g_{\theta}(z)$ represent neural networks. The VAE
architecture is trained via gradient descent by minimizing the evidence lower
bound (ELBO) for each sample $x$.

$log\ p_{\theta}(x)\geq ELBO(\phi,\theta,x)$

$=\mathbf{E}_{z\sim q_{\phi}(z|x)}log\
p_{\theta}(x|z)-KL(q_{\phi}(z|x)||p(z))$

Abstractly, the elements that compose a typical VAE are illustrated at the top
of Figure 1. The term $\mathbf{E}_{z\sim q_{\phi}(z|x)}logp_{\theta}(x|z)$ is
called the reconstruction error or reconstruction loss ($L_{rec}$) of each
sample, often calculated in practice using the mean squared error (MSE)
between the data sample and it’s reconstruction. For a VAE trained on $X\sim
p_{data}$, $L_{rec}$ would be high if the sample is not drawn from $p_{data}$
and would be low for a sample drawn from $p_{data}$. Post training, $L_{rec}$
for any new data sample can be used as an indicator of how well the new sample
conforms to the original data distribution from which the VAE was trained.
Thus, VAEs have been used for anomaly detection tasks in recent years [21, 2,
3].

Figure 1 also illustrates the major processing layers in our implementation of
an eBPF-based container forensic publishing service that employs VAE-based
stability assessment.

Forensics Generation:  The primary source of process forensics are captured
via eBPF and Sysdig, configured to generate forensic records for a selected
set of system calls. Our configuration tracks a default set of 72 Linux system
calls, which represent the most commonly audited system calls [10] [13] to
drive security incident analysis, Host IDS analytics, and application fault
analysis. These systems calls can be categorized into 10 classes of activity:
process events, the set-user-id family of events, network events, file and
directory access events, kernel module load events, process and application
virtualization event, file descriptor replication, file attribute events,
filesystem mount events, and IOCTL events. In addition, the forensic stream
can be augmented with sensor event streams, application log events, system
events, and system resource metrics.

Activity Vector Summarization:  At regular temporal intervals, a forensic
summarization service, or summarizer, computes a vector of statistics from an
analysis of the forensic stream produced during the interval. These statistics
include continuous and categorical measures of events (or even N-grams),
invocation counts, response codes, and other arguments extracted from the
forensic data produced during the interval. The temporal window is
configurable, with a default value of 30 seconds. The activity vector provides
a concise salient behavioral summarization of the runtime pattern observed,
and used by the VAE to assess each runtime interval against its trained data
model (derived from a set of activity vectors produced during the
application’s training period).

Adaptive Forensic Publisher:  We introduce an adaptive forensic publisher that
alters the depth of forensics published based on the VAE’s analysis of the
container activity vector. Here, the depth of forensics published to a remote
forensic repository is selected based on whether the container’s latest
runtime interval is found to be consistent with the vetted-stable forensic
corpus from which the VAE was trained. During intervals in which the VAE-
processed forensic stream produced a low reconstruction error, the VAE latent
model M is published to Elasticsearch. Otherwise, both M and the cached
forensic stream for this interval are published to Elasticsearch. The cache,
depending on the local host configuration, can host forensics across multiple
prior intervals, enabling the publisher to serve Elasticsearch with prior
intervals upon request.

VAE Stability Evaluation: With each incoming summarizer message $x_{in}$ per
container, the publisher checks whether trained normalization and VAE models
for that container exist. They are used to calculate
$L_{rec}=MSE(x_{in},D_{\theta}(E_{\phi}(x_{in})))$, which is then compared to
a pre-set reconstruction error threshold $r_{th}$. $L_{rec}>r_{th}$ is
considered to be a drift in stability, otherwise the operating mode is
considered stable.

VAE Training:  In case a valid container VAE model does not exist, each
incoming vectorized summarizer message is stored until a pre-decided number of
total summarizer messages have been accumulated into a dataset. The dataset
features are normalized and the VAE is trained on the normalized dataset using
mini-batch gradient descent. Post accumulation, both the normalization model
and the VAE neural network model are stored to be used for evaluation.
Important parameters that govern training and evaluation are explained in the
Appendix.

##  III Evaluation

###  III-A Evaluating Multi-container Workload Overhead

We conducted a series of comparative assessments of VAE-based publishing
overhead to that of a standard eBPF-based container monitoring system. The
standard system performs interval caching, and utilizes the Elasticsearch
bulk-API transfer protocol per interval, offering a reasonable and
performance-aware implementation for typical process forensic publishing. We
consider CPU overheard, processing time, network transport overhead, and
Elasticsearch data-storage costs for comparison.

We used the Falco system-call test-generation container provided by the Falco
project [18]. The experiment employed the use of eBPF and Sysdig, configured
to capture 72 Linux system-calls that are commonly used for forensic security
analysis and error diagnosis. The experiment was hosted on an Intel(R)
Xeon(TM) CPU 3.20GHz, 8 core, 16GB DDR2 DRAM, 2TB SATA II Drive, and Docker
hosted on Ubuntu 16. We conducted the performance tests using 10, 20, 30, 40,
and 50 container instances on our host, with 30-second caching intervals, and
then averaged to performance statistics over 50 intervals.

Table 1. Per interval CPU cycles and processing time for all forensic
production and publication computation activity

Table 2. Network transport and Elasticsearch data-storage costs for forensic
production and publication

Table 1 shows the per-interval CPU costs and processing time incurred by the
standard publishing method versus the VAE, for each of five test processing
loads. Perf(1) (Linux command) was used to compute CPU cycles per the two
forensic process algorithms, excluding consideration of the eBPF and sysdig
costs, as these two modules are used and configured identically for both
publishing systems. Of particular interest is the observation that at 50
container instances (90K average syscalls per interval), the standard
publishing method requires approximately 28 seconds to publish the forensics
during the average 30 second interval. In general, when a publisher’s
processing time exceeds the collection interval, this is considered to point
of saturation, as it takes longer to publish the records than to produce them.
In contrast, under the same 50 container workload, the VAE requires only 510
milliseconds to conduct its stability analysis and publish results.

Table 2 shows the network transport overhead and storage costs for the
standard publisher versus the VAE. These costs reflect to overhead involved in
transferring forensic and container-stability assessments to a forensics data
center for incident and fault analysis. The network transport cost is captured
in bytes observed from tcpdump captured from the publisher’s communications
with Elasticsearch over each 5-minute experiment. The table also shows the
average byte storage increase incurred by Elasticsearch during each 5-minute
experiment.

###  III-B CPUMiner: Detecting Container Hijacking using Microservice
Stability Assessment

In addition to performance analyses, we conducted stability-veracity
assessments to assess whether activity vectors captured during stable
operations correspond to a low reconstruction error from a trained model.
Periods in which non-malicious and non-error related operations produce
reconstruction errors above threshold represent incompleteness in the training
set and, depending on their frequency, may indicate a need for enhanced model
retraining. Another assessment perspective is the need to validate that
encounters with operations not represented in the trained model will be
recognized as unstable intervals by the VAE. Our assessment is applicable to a
range of test scenarios, including application hijacking, hardware failures,
infrastructure faults that impact the container, network attacks, and other
fault and attack related scenarios.

Table 3 illustrates an assessment test in which an Nginx-trained VAE observes
a series of runtime intervals during which a malicious CPUMiner[17]
application is integrated into the Nginx container. The assessment begins with
intervals of normal Nginx operation (a shell script that fetches the index
page once per second). As the assessment proceeds, the container experiences a
login shell, followed by the issuing of various shell commands. As these
actions unfold, the VAE’s reconstruction error grows dramatically, by orders
of magnitude as shown in the table. Next, the VAE observes the downloading and
compilation of CPUMiner, a common application that is installed into hijacked
containers to mine Bitcoins. The VAE produces an astronomical reconstruction
error from the initial baseline Nginx operations. Finally, CPUMiner is
initiated in the container, continuing the high reconstruction pattern
indefinitely, or until CPUMiner no longer pollutes the container’s interval
activity vectors.

Activity | VAE reconstruction error
---|---
Normal NGINX processing | 0.014993
Connect with bash shell | 6.298320
Execute shell commands | 19390605381.8
CPUminer package download | 1.3524e+13
Compile CPUminer | 13642381656064
CPUminer execution | 344375512

Table 3. Interval reconstruction errors from a container operating an Nginx-
trained VAE to a gradual progression of actions that are used to download
CPUMiner and install it within the container

###  III-C HTTP Flooding Attack: Detecting Flood Attacks using Microservice
Stability Assessment

Fig 3. Reconstruction errors during a cyclical HTTP-flood attack on the
AtSea_app container

We extend the stability-veracity assessment to a Docker-based demo application
experiencing an HTTP flood attack. We used AtSea [1], a containerized Java
REST application with a database for product inventory, customer data, and
orders, a shopping cart, a Nginx reverse proxy implementing HTTPS and a
payment gateway to simulate certificate management. We used Golang-httpflood
[9] to simulate large number of HTTP requests to overload the AtSea
containers. An instance of the VAE is trained for the atsea_app container and
the Golang-httpflood is invoked cyclically for 1 minute intervals with and 4
minutes of inactivity. Figure 3 shows trained VAE reconstruction errors for
such an attack. The spikes correspond to attack durations.

##  IV Related Work

The topic of anomaly detection in microservices has received recent research
attention. Ohlsson uses Robust Principal Component Analysis (RPCA) techniques
to detect microservice anomalies [15]. Samir et al., describe the use of
Markov Models to detect anomalous container behavior [20] Naikade uses LSTMs
to detect performance anomalies in microservice deployments [14]. Meng et al.
propose detecting variance in application call trees as a strategy to detect
microservice performance anomalies [12]. Magableh uses deep Q-learning
networks (DQNs) to build self-adaptive agents for performance optimization
[11].

Other recent works investigate the applicability of deep learning techniques,
specifically VAEs to the networking domain. Gan et al. [6] and Pol et al. [16]
apply a version of VAEs called conditional VAEs for root cause analysis in
cloud microservices. Xu et al. employ VAEs for anomaly detection in seasonal
KPIs of web applications [21]. For network intrusion detection tasks, a case
study [5] used vanilla auto encoders, while Zavrak et al. used VAEs [22]
trained on flow-based features. Akrami et al. propose beta-VAEs [2] and apply
them to network data for anomaly detection. Other papers have also employed
VAE reconstruction probability as the metric for anomaly detection instead of
the reconstruction error [3]. To our knowledge, our work is the first to
explore the application of VAEs to containerized applications for real-time
stability analysis (as opposed to carrying out anomaly detection on benchmark
networking datasets, as seen in a majority of the aforementioned works) in
conjunction with a dynamic forensic publisher.

##  V Conclusion and Future Work

We present a method to transform container process-level forensics into
interval-based abstract activity vectors, and then input these vectors to a
VAE that produces a metric of container runtime stability. This metric is then
used to drive an adaptive forensic publisher that alters the depth of
information it stores. The approach enables performant determination of the
runtime stability of large numbers of containers and isolation of forensic
data for deeper analysis. A comparative analysis of a VAE-optimized publisher
versus standard process-level forensic auditing shows that our approach
reduces CPU, network, and data storage overhead by orders of magnitude. The
results suggest that fine-grained process-level monitoring can be successfully
scaled to large container ecosystems that are today outside the resource-
consumption costs of traditional forensic data collection services. Future
work includes setting the process-event capturing time window dynamically to
make it difficult for an adversary to spread the attack over multiple windows;
and extending the list of process level events currently being captured for
summarization (details provided in the Appendix).

Appendix

VAE Configuration: The encoder and decoder modules of the VAE are instantiated
each with three hidden layers, each layer has 16 filters in a fully connected
fashion. The dimension of the latent model is chosen to be 10. The network is
trained with the Adam optimizer with learning rate $10^{-4}$, $\beta 1=0.9$,
and $\beta 2=0.999$. The network is trained for 100 epochs. These
hyperparameters were chosen after a grid search on the parameter space for the
specific summarizer dataset per-container. The summarizer dataset for each
container is normalized using a min-max scaler that is also stored with the
VAE model for the purpose of data normalization at test time.

Evaluation Scheme:  We used the perf tool for CPU profiling. Perf stats are
captured by attaching the tool to the process IDs of Sysdig chisel,
summarizer, and publisher processes running on the Linux test machine.
Networking specifics are extracted using a running tcpdump process. Storage
usage metrics are extracted from elastic search running on a separate server.

Reconstruction error and reconstruction threshold: We employ VAE
reconstruction error for each new sample as an indicator of stability.
Throughout training, the VAEs minimize the reconstruction error on a dataset
as a means of learning the data distribution. A sample training curve for a
VAE trained on a dataset is represented in Fig. 4.

Fig 4. Reconstruction errors during training. VAE train reconstruction error
settles to $r_{last}$ after epoch 40

The training reconstruction error settles to a minimum value $r_{last}$ after
training. For each new sample, deviation from $r_{last}$ signifies that the
data sample is different from the original training data or is out of
distribution. Owing to this fact, we set the reconstruction threshold $r_{th}$
in two major ways for our experiments:

  * •

Based on a heuristic/experience: This approach sets $r_{th}$ based on observed
reconstruction errors generated during an instable period.

  * •

$k$ standard deviations away from the mean : This approach records the mean
($r_{mean}$) and standard deviation ($r_{sd}$) of training reconstruction
errors in the last training iteration. The reconstruction threshold is then
set using the following formula : $r_{th}=r_{mean}+k*r_{sd}$ where $k$ is set
by the user. This approach has the benefit of making the VAE more flexible or
conservative based on the value of $k$.

Summarization of Forensic Events:  We used a (non-exhaustive) list of the
following kernel events:

  * •

Process Events:

    * –

common fork __x64_sys_fork/ptregs

    * –

common vfork __x64_sys_vfork/ptregs

    * –

64 execve __x64_sys_execve/ptregs

    * –

64 execveat __x64_sys_execveat/ptregs

    * –

common exit __x64_sys_exit

    * –

common kill __x64_sys_kill

    * –

64 ptrace __x64_sys_ptrace

    * –

common prctl __x64_sys_prctl

    * –

common arch_prctl __x64_sys_arch_prctl

    * –

x32 execve __x32_compat_sys_execve/ptregs

    * –

x32 ptrace __x32_compat_sys_execveat/ptregs

  * •

Set User ID Family Events:

    * –

common setuid __x64_sys_setuid

    * –

common setgid __x64_sys_setgid

    * –

common setpgid __x64_sys_setpgid

    * –

common setsid __x64_sys_setsid

    * –

common setreuid __x64_sys_setreuid

    * –

common setregid __x64_sys_setregid

    * –

common setgroups __x64_sys_setgroups

    * –

common setresuid __x64_sys_setresuid

    * –

common setresgid __x64_sys_setresgid

    * –

common setfsuid __x64_sys_setfsuid

    * –

common setfsgid __x64_sys_setfsgid

  * •

Network Events:

    * –

common socket __x64_sys_socket

    * –

common connect __x64_sys_connect

    * –

common accept __x64_sys_accept

    * –

common bind __x64_sys_bind

    * –

common listen __x64_sys_listen

    * –

common accept4 __x64_sys_accept4

  * •

File and Directory Access Events:

    * –

common open __x64_sys_open

    * –

common close __x64_sys_close

    * –

common openat __x64_sys_openat

    * –

common mkdir __x64_sys_mkdir

    * –

common rmdir __x64_sys_rmdir

    * –

common rename __x64_sys_rename

    * –

common creat __x64_sys_creat

    * –

common link __x64_sys_link

    * –

common unlink __x64_sys_unlink

    * –

common symlink __x64_sys_symlink

  * •

File and Directory Access Events (cont.)

    * –

common mknod __x64_sys_mknod

    * –

common mkdirat __x64_sys_mkdirat

    * –

common mknodat __x64_sys_mknodat

    * –

common unlinkat __x64_sys_unlinkat

    * –

common renameat __x64_sys_renameat

    * –

common linkat __x64_sys_linkat

    * –

common symlinkat __x64_sys_symlinkat

    * –

common fchmodat __x64_sys_fchmodat

    * –

common renameat2 __x64_sys_renameat2

  * •

Kernel Module Load Events:

    * –

64 create_module

    * –

common init_module __x64_sys_init_module

    * –

common delete_module __x64_sys_delete_module

    * –

x32 kexec_load __x32_compat_sys_kexec_load

  * •

Process/App Virtualization Event:

    * –

common clone __x64_sys_clone/ptregs

    * –

common clone3 __x64_sys_clone3/ptregs

    * –

process_vm_readv

      * *

64 __x64_sys_process_vm_readv

      * *

x32 __x32_compat_sys_process_vm_readv

    * –

process_vm_writev

      * *

64 __x64_sys_process_vm_writev

      * *

x32 __x32_compat_sys_process_vm_writev

  * •

File Descriptor Replication:

    * –

common dup __x64_sys_dup

    * –

common dup2 __x64_sys_dup2

    * –

common dup3 __x64_sys_dup3

  * •

File Attribute Events:

    * –

common chmod __x64_sys_chmod

    * –

common fchmod __x64_sys_fchmod

    * –

common chown __x64_sys_chown

    * –

common fchown __x64_sys_fchown

    * –

common lchown __x64_sys_lchown

    * –

common fchownat __x64_sys_fchownat

  * •

Filesystem Mount Events:

    * –

common mount __x64_sys_mount

    * –

common umount2 __x64_sys_umount

    * –

common fsmount __x64_sys_fsmount

  * •

IOCTL Events:

    * –

64 ioctl __x64_sys_ioctl

    * –

x32 ioctl __x32_compat_sys_ioctl

These event indicators are extracted from the eBPF forensic log and passed
through our summarizer module to be distilled in a numerical vector format. A
collection of these vectors form our dataset.

## References

  * [1] A. Akkuş, “The atsea sample shop,” https://github.com/dockersamples/atsea-sample-shop-app.
  * [2] H. Akrami, A. Joshi, and S. Aydore, “Robust variational autoencoder for tabular data with beta divergence,” _ICML Workshop_ , 2020.
  * [3] J. An and S. Cho, “Variational autoencoder based anomaly detection using reconstruction probability,” _Thesis_ , 2014.
  * [4] F. G. et.al, “IDC Futurescape: Worldwide IT industry 2020 predictions,” _IDC Research_ , vol. US45599219, October 2019, https://www.idc.com/getdoc.jsp?containerId=US45599219.
  * [5] G. Fernandez and S. Xu, “A case study on using deep learning for network intrusion detection,” _arxiv_ , 2019.
  * [6] Y. Gan, D. Lo, S. Dev, and C. Delimitrou, “Sage: Leveraging ml to diagnose unpredictable performance in cloud microservices,” _Machine Learning for computer architecture and systems_ , 2020.
  * [7] S. Incorporated, “Welcome to Sysdig!” https://github.com/draios/sysdig.
  * [8] D. Kingma and M. Welling, “Auto encoding variational bayes,” _International Conference on Learning Representations_ , 2014.
  * [9] Leeon123, “The golang httpflood,” https://github.com/Leeon123/golang-httpflood.
  * [10] U. Lindqvist and P. A. Porras, “expert-bsm: A host-based intrusion detection solution for sun solaris,” in _17th Annual Computer Security Applications Conference (ACSAC 2001), 11-14 December 2001, New Orleans, Louisiana, USA_. IEEE Computer Society, 2001.
  * [11] B. Magableh, “Deep q learning for self-adaptive distributed microservices architecture,” _IEEE Acess_ , 2019.
  * [12] L. Meng, F. Ji, Y. Sun, and T. Wang, “Detecting anomalies in microservices with execution trace comparison,” _Future Generation Computer Systems_ , vol. 116, pp. 291–301, 2021.
  * [13] A. Mounji and B. L. Charlier, “Tools for intrusion detection: Results and lessons learned from the asax,” in _The First Recent Advances in Intrusion Detection, International Workshop, RAID 1998, Louvain-la-Neuve, Belgium), year = 1998,_.
  * [14] P. P. Naikade, “Automated anomaly detection and localization system for a microservices-based cloud system,” _Electronic Thesis and Dissertation Repository_ , 2020.
  * [15] J. OHLSSON, “Anomaly detection in microservice infrastructures,” _M.S. Thesis_ , 2018.
  * [16] A. Pol, V. Berger, G. Cerminara, C. Germain, and M. Pierini, “Anomaly detection with conditional variational autoencoders,” _ICMLA_ , 2019.
  * [17] Pooler, “CPUMiner: a cpu miner for litecoin, bitcoin, and other cryptocurrencies,” https://sourceforge.net/projects/cpuminer/.
  * [18] T. F. Project, “Falco: Cloud native runtime security,” https://falco.org/.
  * [19] L. Pustina, “Sysdig Chisel,” https://github.com/lukaspustina/sysdig_chisels.
  * [20] A. Samir and C. Pahl, “Dla: Detecting and localizing anomalies in containerized microservice architectures using markov models,” in _2019 7th International Conference on Future Internet of Things and Cloud (FiCloud)_ , 2019, pp. 205–213.
  * [21] H. Xu, W. Chen, N. Zhao, Z. Li, J. Bu, Z. Li, Y. Liu, Y. Zhao, D. Pei, Y. Fend, J. Chen, and H. Wang, Z mand Qiao, “Unsupervised anomaly detection via variational auto-encoder for seasonal kpis in web applications,” _World Wide Web Conference_ , 2018.
  * [22] S. Zavrak and M. İskefiyeli, “Anomaly-based intrusion detection from network flow features using variational autoencoder,” _IEEE Access_ , 2020.