[orcid=0000-0003-1874-1854] [1]

[orcid=0000-0001-5132-9488]

[orcid=0000-0003-1129-5038]

[orcid=0000-0002-1295-7657]

1]organization=Department of Engineering and Applied Science (DISA),
University of Bergamo, addressline=Via Pasubio 3, city=Dalmine (BG),
postcode=24044, state=Italy,

2]organization=Department of Economics, University of Bergamo, addressline=Via
dei Caniana 2, city=Bergamo, postcode=24127, state=Italy,

[1]Corresponding author

# Scenario analysis of livestock-related PM2.5 pollution based on
heteroskedastic geostatistical modelling

Jacopo Rodeschini<EMAIL_ADDRESS>Alessandro Fassò
<EMAIL_ADDRESS>Alessandro Fusta Moro<EMAIL_ADDRESS>Francesco Finazzi<EMAIL_ADDRESS>[  [

###### Abstract

The air in the Lombardy region, Italy, is one of the most polluted in Europe
because of limited air circulation and high emissions levels. There is a large
scientific consensus that the agricultural sector has a major impact on air
quality. In Lombardy, livestock activities are widely acknowledged to be
responsible for approximately 97% of regional ammonia emissions due to the
high density of livestock. The main objective of our study is to quantify the
relationship between ammonia emissions and PM2.5 concentrations in the
Lombardy region and evaluate PM2.5 changes due to the reduction of ammonia
emissions through scenario analysis. In particular, the study refers to the
years between 2016 and 2020 inclusive. The information contained in the data
is exploited using a spatiotemporal model capable of handling spatial and
temporal correlation, as well as missing data. In this study, we propose a
heteroskedastic extension of the Hidden Dynamic Geostatistical Model (HDGM)
which is a two-level hierarchical model suitable for complex environmental
processes. Scenario analysis will be carried out on high-resolution maps of
the Lombardy region showing the changes in PM2.5 across the area. As a result,
it is shown that a 26% reduction in $NH_{3}$ emissions in the wintertime could
reduce the PM2.5 average by $2.09$ $\mu g/m^{3}$ while a 50% reduction could
reduce the PM2.5 average by $4.02$ $\mu g/m^{3}$ which correspond to a
reduction close to 5% and 10% respectively. Finally, results are detailed by
province and land type.

###### keywords:

PM2.5 concentrations Scenario analysis Ammonia emissions Air quality Lombardy
Italy

##  1 Introduction

Air pollutants can be divided into two categories: primary and secondary.
Primary pollutants are those that are directly released into the atmosphere,
while secondary pollutants are formed in the atmosphere from precursor gases
through chemical reactions and microphysical processes. Ammonia (NH3) is a key
precursor gas for secondary particulate matter (PM) such as PM10 and PM2.5. In
the Lombardy region, livestock and fertilisers are the main sources of
ammonia, close to 97% of the overall regional emissions [14, 13, 3]. The total
annual NH3 emissions are roughly constant during the considered period [14].
However, seasonal practises of agricultural activities determine strong
variations within the year [13]. In the Po basin, PM2.5 composition in winter
has a very significant secondary component, usually exceeding 50% with the
relevant part composed of ammonium salts resulting from chemical reactions
between ammonia and other pollutants. This process is favoured during
accumulation episodes under unfavourable meteorological conditions for
dispersion. Indeed, due to its stability in air and limited air circulation in
the Po Valley, PM2.5 accumulates often at hazardous concentrations and has a
negative impact on human health. This is especially true during the cold
months when a reduced mixing boundary layer is combined with intense emissions
from endothermic-powered vehicles and domestic heating systems. As a result of
this high level of pollution, Lombardy has been found to have one of the
highest PM2.5-related mortality rates among all regions of Europe, with
100-150 premature deaths per 100.000 inhabitants [8].

Several studies [16, 6, 19] found that in Europe the reduction of $NH_{3}$
emissions is one of the most effective control strategies used to mitigate
PM2.5 in both summer and winter, mainly due to a significant decrease in
ammonium salts. In [16] five control strategies were examined separately, and
a 50 % reduction of gaseous emissions (SO2, NH3, NOx, anthropogenic VOCs) are
tested to explore the sensitivity of concentrations to emissions. They showed
that in most European areas the reduction of NH3 emissions during winter and
summer periods is more effective in reducing total PM2.5 than reductions of
other gas precursors. In [6] a 50% of ammonia emissions reduction was applied
specifically for the Lombardy region resulting in a reduction of total PM2.5
up to 2.4 $\mu g/m^{3}$. Furthermore, the efficiency of PM2.5 controls
strategy is affected by the nonlinearity behaviour of the sulphate-nitrate-
ammonia system, as shown by [19].

In Italy, legal levels for air pollutants are regulated by Legislative Decree
No. 155 of Aug. 13, 2010 (D.Lgs. 155 /2010), which implements European
Directive 2008 /50 /EC on air quality and cleaner air for Europe. In
particular, D.Lgs. 155 /2010 assigns the regions the task of adopting a plan
with the necessary measures to act on the main emissions sources to achieve
the limit values. Regarding $PM_{2.5}$, the limit value is identified as the
annual average below 25 $\mu g/m^{3}$, far from the value of 5 $\mu g/m^{3}$
identified by the last report of WHO as a safe threshold. Many cities in
northern Italy have not complied with national limits, leading the European
Commission to initiate infringement proceedings (after those on ozone and
$PM_{10}$) and forcing Italian regions to increase efforts to improve air
quality. In Lombardy, the Regional Plan for Air Quality Interventions (PRIA)
envisages the achievement of air quality objectives in 2025, and in the latest
update (2021) the reductions are in line with expectations.

This work aims to use scenario analysis to assess the reductions in PM2.5
concentration due to ammonia reductions. A geostatistical model is used to
account for the spatiotemporal correlation. In the last decades, the
development of spatiotemporal statistical models of air quality has been
extensive (see, e.g., [2, 5, 20]). Chemical transport model (CTM) simulations
have been applied by [22, 3, 19] to study ammonia emissions reduction
scenarios. In [21] an ensemble strategy based on three CTM is adopted to model
the AQ in Lombardy. An alternative approach is developed in [9] where they
study the sensitivity of PM2.5 to $NH_{3}$ concentrations. Focusing on the
Lombardy region, they consider a seasonal linear model with serial
correlation, correlation among stations, and heteroskedasticity. In this
study, we propose a heteroskedastic extension of the Hidden Dynamic
Geostatistical Model (HDGM) [2] which is a two-level hierarchical model
suitable for complex environmental processes [24]. HDGM has been compared with
the Generalised Additive Mixed Model (GAMM) and a Random Forest model Residual
Kriging (RFRK) in PM2.5 modelling in Lombardy and has demonstrated better
performance in cross-validation [18].

The rest of the paper is organised as follows. Section 2 discusses the
methodology adopted. In particular Section 2.1 introduces the heteroskedastic
extension of the HDGM model. We propose a two-step procedure to estimate the
parameters. In the first step, the maximum likelihood (MLE) method is used to
estimate the deterministic trend and the autocorrelation structure. In the
second step, the error variances are estimated on the first-step residuals.
Section 2.2 defines the scenario analysis approach. Section 3 discusses the
application of the HDGM model to the observed daily AQ data for the Lombardy
Region, between 2016 and 2020. More precisely, we consider ground PM2.5 as the
response variable, while weather, emissions flows, and land use are considered
explanatory variables. Section 3.1 describes the dataset and the data
explored. Section 3.2 discussed the specification of the models while Section
3.2 reported the residual diagnostics. Finally, a scenario analysis of Section
3.4 is performed to assess the PM2.5 changes due to NH3 emissions on high-
resolution maps. The Conclusion closes the paper.

##  2 Methodology

In this section, we detail the methods adopted to perform scenario analysis.
To perform a scenario analysis, it is important to predict the PM2.5
concentrations over the Lombardy surface. Although the HDGM model is able to
handle multivariate spatiotemporal data, in this study we refer to the
univariate case. The next section summarises a statistical modelling approach
to estimate the PM2.5 concentration starting from the data observed in the AQ
monitoring network.

###  2.1 Heteroskedastic Hidden Dynamic Geostatistical Model

To understand the relationship between emissions and particulate matter,
taking into account spatial and temporal correlation, we propose a
heteroskedastic extension of the Hidden Dynamic Geostatistical Model (HDGM)
[4, 2] which is a two-level hierarchical model. The hierarchy is constructed
by putting together two conditional submodels. At the first level, the
observation variability is modelled by the measurement equation, which is
essentially given by a deterministic trend, a stochastic signal, and an error.
The signal is defined at the second level of the hierarchy and handles the
spatiotemporal correlation through a Markovian process and an innovation term
which is a zero-mean Gaussian process with a spatial exponential covariance
function. A two-step procedure estimates the parameters. In the first step,
the maximum likelihood (MLE) method estimates the deterministic trend and the
autocorrelation structure. In the second step, the error variances are
estimated on the basis of the first-step residuals.

Let $y(\bm{s},t)$ be the response variable observed at site
$\bm{s}\in\mathcal{S}$ and discrete time $t=1,...,T$ where $\mathcal{S}$ is
the geographic region of interest. The heteroskedastic HDGM is defined as
follows:

$\begin{cases}y(\bm{s},t)=\bm{x}(\bm{s},t)\bm{\beta}+vz(\bm{s},t)+\epsilon(\bm{s},t)\\\
z(\bm{s},t)=gz(\bm{s},t-1)+\eta(\bm{s},t)\end{cases}$ (1)

which under general assumptions can be considered as a classical state-space
model [7]. The $\bm{\beta}$ is a vector of fixed effects; $\bm{x}(\bm{s},t)$
is a vector of covariates that accounts for all exogenous effects; the
measurement error $\epsilon(\bm{s},t)\sim N(0,\sigma^{2}_{\epsilon}(t))$ which
independent in space and time; $z(\bm{s},t)$ is a Markovian univariate process
ruled by the transition coefficients $g$; $v$ is a scaling coefficient; the
innovation term $\eta(\bm{s},t)$ is a zero-mean Gaussian process
$GP(0,\rho(\|\bm{s}-\bm{s}^{\prime}\|;\theta))$ independent in time where
$\rho$ is the exponential correlation function and
$\|\bm{s}-\bm{s}^{\prime}\|$ is the distance between $\bm{s}$ and
$\bm{s}^{\prime}\in\mathcal{S}$. The model parameters set is
$\bm{\Phi}=\\{\bm{\beta},g,\theta,v,\sigma^{2}_{\epsilon}(t)\\}$.

A two-step procedure is performed to estimate the parameters:

  1. 1.

First-step: the parameters are estimated by the maximum likelihood (MLE)
method assuming errors with constant variance, or equivalently
$\epsilon(\bm{s},t)\sim N(0,\sigma^{2}_{\epsilon})$ $IID$. To do this, we use
the open-source software D-STEM [23].

  2. 2.

Second-step: empirical estimation of the time-dependent error variance, say
$\hat{\sigma}^{2}_{\epsilon}(t)$, based on the first-step residuals.

The first-step residuals $e(\bm{s},t)$ are computed as:

$e(\bm{s},t)=y(\bm{s},t)-\hat{y}(\bm{s},t)$ (2)

where $y(\bm{s},t)$ and $\hat{y}(\bm{s},t)$ are the observation and prediction
in the location $\bm{s}$ at time $t$ considering the MLE of the parameters in
the first step. Prediction $\hat{y}(\bm{s},t)$ at site $\bm{s}$ and time
$t=1,...,T$ is calculated using the plug-in approach as follows:

$\hat{y}(\bm{s},t)=\bm{x}(\bm{s},t)\hat{\bm{\beta}}+\bm{z_{t}}^{T}(\bm{s})$
(3)

where $\bm{z_{t}}^{T}(\bm{s})$ is the estimate of the latent variable done by
the Kalman smoother algorithm.

In the second step, since we expect a high daily variability, we estimate
$\hat{\sigma}_{\epsilon}(t)$ using a smoothing function of the daily
$RMSE_{t}$. More precisely, assuming the residual mean equal to zero,
$\hat{\sigma}_{\epsilon}(t)=f(RMSE_{t})$ where $f$ is an interpolation spline
over the monthly averages of $RMSE_{t}$ where the monthly averages are
associated with the first day of the month. The daily $RMSE_{t}$ is computed
as:

$RMSE_{t}=\sqrt{\frac{\sum_{\bm{s}\in\mathcal{S}}e(\bm{s},t)^{2}}{|\mathcal{S}|}}$
(4)

where $|\mathcal{S}|$ denote the number of sites involved in the computation

###  2.2 Scenario analysis approach

This section presents the scenario analysis approach adopted to assess the
PM2.5 sensitivity to the NH3 emission. We consider the "What-if" scenario
analysis approach: "What would happen if the NH3 reductions were fully
implemented at time 0 (1st January 2016)". We compare the baseline case of
observed NH3 emissions, denoted as $\hat{y}(\bm{s},t)$, with a prediction in
which the NH3 emissions are reduced by a $\alpha$ factor,
$\hat{y}^{\alpha}(\bm{s},t)$. The daily PM2.5 changes $\Delta_{y}(\bm{s},t)$
and the associated uncertainty can be calculated as:

$\begin{cases}\Delta_{y}(\bm{s},t)=\hat{y}(\bm{s},t)-\hat{y}^{\alpha}(\bm{s},t)=\Delta_{\bm{x}}\hat{\bm{\beta}}\\\
MSE[\Delta_{y}(\bm{s},t)]=\Delta_{\bm{x}}^{\prime}\bm{\Sigma}_{\hat{\bm{\beta}}}\Delta_{\bm{x}}+V[\epsilon_{t}-\epsilon_{t}^{\alpha}]\end{cases}$
(5)

where $\Delta_{\bm{x}}=\bm{x}(\bm{s},t)-\bm{x}^{\alpha}(\bm{s},t)$;
$\hat{\bm{\beta}}$ and $\bm{\Sigma}_{\hat{\bm{\beta}}}$ are the estimates of
fixed effects and the associate covariance matrix common for both predictions;
$\epsilon_{t}$ and $\epsilon_{t}^{\alpha}$ are the measurement errors of both
predictions. Once we predict the daily $\Delta_{y}(\bm{s},t)$ on a regular
grid, we compute the average of the daily $\Delta_{y}(\bm{s},t)$ over space
subset $\mathcal{D}$ and time interval $\mathcal{I}$. Let $|D|$ and $|I|$ be
the number of pixels and days involved in the average computation. Assuming
that $\epsilon_{t}$ and $\epsilon_{t}^{\alpha}$ are independent with the same
variance, we compute the daily average
$\bar{\Delta}_{y}(\mathcal{D},\mathcal{I})$ as:

$\begin{cases}\bar{\Delta}_{y}=\frac{1}{|D||I|}\sum_{\bm{s}\in\mathcal{D}}\sum_{t\in\mathcal{I}}\Delta_{y}(\bm{s},t)=\bar{\Delta}_{\bm{x}}\hat{\bm{\beta}}\\\
MSE[\bar{\Delta}_{y}]\leq\bar{\Delta}_{\bm{x}}^{\prime}\bm{\Sigma}_{\hat{\bm{\beta}}}\bar{\Delta}_{\bm{x}}+\frac{2}{|D|^{2}|I|^{2}}\sum_{\bm{s}\in\mathcal{D}}\sum_{t\in\mathcal{I}}\hat{\sigma}_{\epsilon}^{2}(s,t)\end{cases}$
(6)

where $\hat{\sigma}_{\epsilon}^{2}(\bm{s},t)$ is the estimate of the variance
of the measurement error and $\bar{\Delta}_{\bm{x}}$ is the average of the
design matrix. Derivation of Eq. 6 is reported in Appendix A.

According to [21] we propose two different scenarios for NH3 emissions
reduction: (i) Scenario named PRIA that accounts for a reduction of 26% in NH3
emission; (ii) Scenario named Strong that is characterised by a strong
reduction of 50% in ammonia emissions. The PRIA scenario is based on the PRIA
plan, which identifies the actions needed to reduce ammonia emissions by 26%.
Table 1 summarises the two proposed scenarios. The average and uncertainty of
the PM2.5 changes are shown on high-resolution maps using a regular grid (0.1∘
$\times$ 0.1∘) over the studied area.

Table 1: Proposed scenarios on ammonia reduction. The scenarios are implemented by setting $\alpha$ of Eq. 5 according to NH3 reduction (%) defined by each scenario. Scenarios | NH3 reduction (%) | $\alpha$
---|---|---
PRIA | 26 | 0.26
Strong | 50 | 0.50

##  3 Lombardy case study

In this section, we use the methodology discussed in Section 2 to assess the
role of the $NH_{3}$ emissions on the $PM_{2.5}$ concentration in the Lombardy
region.

###  3.1 Data description

Our research is based on the Agrimonia dataset [11], which is a comprehensive
daily spatiotemporal dataset for AQ modelling of the Lombardy region. This
dataset is open access and is available through the Zenodo repository [10].
Covers the period 2016 to 2021 and includes daily concentrations of air
pollutants, weather conditions, emissions flows, land use, and livestock
densities. Since geostatistical methods can benefit from neighbouring
information to improve performances close to the borders, the dataset provides
data also for an area around the Lombardy region, obtained by applying a 0.3∘
buffer around the regional borders. From the dataset, we consider the period
between 2016 to 2020 and the following variables: PM2.5 measured at 45 ground
stations (30 of them belonging to the Lombardy area, while the remaining 15
are located in the neighbouring area), wind speed (average wind speed at 100
$m$), temperature (air temperature at 2 $m$), relative humidity (RH), the
height of the boundary layer (BLH, daily maximum height of the air mixing
layer), high vegetation index (HVI, high vegetation abundance), sulphur
dioxides emissions ($SO_{2}$), agriculture ammonia emissions ($NH_{3}$, sum of
ammonia emissions from agricultural waste burning, livestock manure management
and agriculture soil) and nitrogen oxide emissions ($NO_{x}$). Table 2
summarises the numerical variables from the Agrimonia dataset used in this
study. Furthermore, we create new binary variables: $Rain$ if total
precipitation exceeds the threshold of 1 $mm$, Urban if the stations belong to
urban land type, and new categorical variables Season and Month.

Table 2: Numerical variables selected from the Agrimonia dataset [11] and main descriptive statistics. Name | min | mean | std | max
---|---|---|---|---
PM2.5 $[\mu g/m^{3}]$ | 1.00 | 21.58 | 16.63 | 182.00
Wind speed $[m/s]$ | 0.56 | 2.59 | 1.33 | 11.93
Temperature $[C^{\circ}]$ | -11.94 | 13.44 | 7.98 | 32.88
RH $[\%]$ | 25.87 | 74.34 | 12.26 | 99.10
BLH $[m]$ | 27.38 | 1043.91 | 556.77 | 4420
HVI $[m^{2}/m^{2}]$ | 0.86 | 2.23 | 0.70 | 5.03
SO2 $[t/km^{2}]$ | 0.05 | 3.34 | 5.53 | 45.06
NH3 $[t/km^{2}]$ | 0.13 | 11.55 | 12 | 72.03
NOx $[t/km^{2}]$ | 0.75 | 14.06 | 16.73 | 118.60

Evidence from [18, 13] suggests that there is a seasonal pattern of higher
concentrations of PM2.5 in winter due to meteorological conditions that reduce
air circulation and create clusters of similar concentrations in neighbouring
areas, implying a strong spatial and temporal dependence. The spatiotemporal
variogram [4] in Figure 1 confirms this, as it shows that the variance
increases with distance in both space and time, with the spatial correlation
stabilising at around 90 $km$ and the temporal correlation at around 7 days.
The variogram lines are roughly parallel, which is consistent with a separable
autocovariance structure. The temporal dynamics are further investigated
through the Partial Auto-Correlation Function (PACF) which highlights the
first lag (1 day) as the most significant auto-correlated component,
indicating that a model with first-order Markovian dynamics is appropriate.
The magnitude at the first lag is $\sim 0.8$ suggesting that PM2.5 is
relatively stable in the atmosphere and has a similar behaviour throughout the
region. Details of PACF are discussed in the Appendix B.

Figure 1: Spatiotemporal variogram [4] computed on PM2.5 daily observations at
45 stations from 2016 to 2020.

A special focus is dedicated to the NH3 emissions. The emissions flow included
in the Agrimonia dataset, harmonised at a daily temporal resolution, [11] were
provided by the Copernicus Atmosphere Monitoring Service (CAMS) [12] on a
monthly basis. The geographical distribution of the annual NH3 emissions in
2020 is shown in Figure 2 while Table 3 summarises the annual $NH_{3}$
emissions by province. In particular, the provinces of Brescia, Cremona, and
Mantova are characterised by higher emissions. It is worth noting, that the
air quality monitoring network is not specifically designed to monitor ammonia
emissions, indeed many monitoring stations are located in areas with low
emissions levels and no one is located where the NH3 emissions reach the major
peaks, as shown in Figure 2.

Figure 2: NH3 emissions density ($t/km^{2}$) over the augmented Lombardy region. The AQ stations (cyan circles) are spread randomly over the studied area. The continuous blue line represents the boundary of Lombardy while the dotted ones represent the provinces.  Table 3: Total NH3 emissions by province in 2020. Province | NH3 $[t]$
---|---
Bergamo (BG) | 7546
Brescia (BS) | 29524
Como (CO) | 1517
Cremona (CR) | 16892
Lecco (LC) | 615
Lodi (LO) | 6552
Mantova (MN) | 20336
Milano (MI) | 14013
Monza e della Brianza (MB) | 714
Pavia (PV) | 8199
Sondrio (SO) | 4601
Varese (VA) | 1052

Another important aspect is the seasonality and the composition regimes in the
atmosphere. Table 4 summarises the daily average of temperature, emissions
flows, and $PM_{2.5}$ by season. In particular, winter is characterised by
higher PM2.5 concentrations. This is related to the lower boundary layer
height in winter, causing less dispersion of airborne pollution and increasing
the importance of precursors.

As the objective of our study is to characterise the role of agricultural
emissions and specifically NH3 emissions on PM2.5 pollution, here we focus
only on the NH3 reduction.

Table 4: Seasonal average of daily Temperature $[C^{\circ}]$, $NH_{3}$, $NO_{x}$ and $SO_{2}$ emissions density $[mg/m^{2}]$ and $PM_{2.5}$ concentrations $[\mu g/m^{3}]$. Season | Temperature | $NH_{3}$ | $NO_{x}$ | $SO_{2}$ | $PM_{2.5}$
---|---|---|---|---|---
Winter | 3.88 | 5.21 | 19.97 | 4.11 | 35.87
Spring | 12.67 | 17.45 | 11.33 | 2.98 | 16.54
Summer | 23.23 | 14.71 | 10.45 | 2.80 | 12.19
Autumn | 13.83 | 8.70 | 14.59 | 3.46 | 21.62

###  3.2 Model estimation and selection

For the purposes of this study, we consider the univariate case of the HDGM
model where the response variable is the PM2.5 in the original scale. We can
use the Wilkinson notation to define the fixed-effects terms as follows:

$PM_{2.5}\sim(Intercept)+Weather+Land\\\ +Emission:Seasonality$ (7)

where Emissions includes NH3, NOx and SO2 and Seasonality one of the two
between Season and Month. To simplify the scenario implementation, we consider
all variables in their original scale and the HDGM coefficients can be
interpreted as marginal (linear) effects. Transformations such as a
logarithmic transformation of the response variable did not lead to better
prediction results and model fits. Although prediction is not the main
objective of this paper, we evaluate the model fit in cross-validation (CV).
In particular, we use the leave-one-station-out cross-validation scheme
(LOSOCV), a variation of the commonly used leave-one-out cross-validation
approach applied in the spatiotemporal framework [17]. This method withholds a
complete time series of a single station from the training of the model and
uses it to evaluate the model’s prediction performance. As the objective of
this study is to characterise the model over the Lombardy region, only the
stations belonging to Lombardy are involved in the CV procedure (30 stations)
while the remaining (15 stations) are always in the training set. After the
implementation of the LOSOCV approach, we obtained prediction results for the
30 stations included in the validation process. To specify the model, we
compare the cross-validation Root Mean Square Error (CV-RMSE) and Akaike
information criterion (AIC) of 32 candidate models, which are characterised by
different subsets of predictors and their interactions. The CV-RMSE is
obtained as an average of the RMSE computed by the station. We specifically
designed the model selection to explore the connection between the source of
emissions and the type of seasonality. To do this, we tested the model with
various interactions between emission and seasonality, as seen in Eq. 7.
According to the equifinality concept [1, 15], we decided to avoid automatic
feature selection. In particular, we identified a set of acceptable models
instead of a single best and focused on nearly the best models, in which the
above scores are very close to the best. Then we chose the model that most
closely addresses our research questions. Figure 3 shows the comparison
between the 32 candidate models with respect to CV-RMSE and AIC sorted
according to the number of covariates. As a result of model selection, we
select the model $\\#$12\. Details of each model are provided in the
Supplementary Material (S1).

Figure 3: CV-RMSE (blue circle) with confident interval (blue bar) and AIC
(magenta cross) comparison of the 32 candidate models sorted on the number of
covariates. The chosen model (#12) is highlighted with a vertical dashed line.

###  3.3 Model diagnostics

The selected model ($\\#$12) includes 17 covariates as follows:

$PM_{2.5}\sim(Intercept)+Windspeed+Temperature\\\
+RH+Rain+BLH+Urban+HVI+SO_{2}\\\ +NO_{x}+NH_{3}+(NO_{x}+NH_{3}):Season$ (8)

The model reaches $AIC=Best+0.23\%$ (where Best identifies the overall
absolute minimum AIC) and $CV-RMSE=Best+2\%$ (where Best identifies the
overall absolute minimum CV-RMSE), resulting in an increasing CV-RMSE of
$0.12\mu g/m^{3}$.

Table 5 summarises the estimated parameters $\bm{\beta}$ of the linear
regression model used for the large-scale component. All coefficients are
statistically significant. The signs of most of the coefficients are
consistent with our expectations. For example, every 1 $m/s$ of extra wind
speed is related to a decrease of 1.67 $\mu g/m^{3}$ of PM2.5. Relative
humidity is positively related, so high humidity levels (100%) are associated
with an increase of 0.14 $\mu g/m^{3}$. Maximum BLH is negatively associated,
each increase of 1000 $m$ is related to a decrease of 3.6 $\mu g/m^{3}$ of
PM2.5. Table 6 shows the average effect on $PM_{2.5}$ due to temperature,
$NH_{3}$, $NO_{x}$ and $SO_{2}$. As expected, the role of $NH_{3}$ is more
relevant in the winter.

Table 5: Fixed effect coefficients $\bm{\hat{\beta}}$ estimates. Name | $\beta$ | std | |t| | p-value
---|---|---|---|---
(Intercept) | 26.34 | 1.24 | 21.18 | 0
Wind speed | -1.67 | 0.05 | 32.25 | 0
Temperature | -0.44 | 0.04 | 10.02 | 0
RH | 0.14 | 0.01 | 15.44 | 0
BLH | -3.60e-03 | 0.00 | 19.86 | 0
Rain | -0.74 | 0.08 | 9.09 | 0
Urban | -1.35 | 0.11 | 12.51 | 0
HVI | -1.24 | 0.14 | 9.02 | 0
SO2 | -0.11 | 0.02 | 5.78 | 7.43e-09
NOx:Winter | 0.09 | 0.01 | 7.53 | 4.94e-14
NOx:Summer | -0.06 | 0.02 | 3.19 | 1.40e-01
NOx:Spring | -0.05 | 0.02 | 2.81 | 0.47
NOx | 0.10 | 0.01 | 9.29 | 0
NH3:Winter | 0.75 | 0.05 | 16.15 | 0
NH3:Summer | -0.10 | 0.03 | 3.84 | 1.24e-04
NH3:Spring | -0.15 | 0.03 | 5.58 | 2.39e-08
NH3 | 0.15 | 0.02 | 6.25 | 4.08e-10
Table 6: Average season effect on $PM_{2.5}$ due to the Temperature, $NH_{3}$, $NO_{x}$ and $SO_{2}$ in $[\mu g/m^{3}]$ Season | Temperature | $NH_{3}$ | $NO_{x}$ | $SO_{2}$ | $PM_{2.5}$
---|---|---|---|---|---
Autumn | -6.08 | 1.31 | 1.46 | -0.38 | 21.62
Spring | -5.57 | 0.00 | 0.57 | -0.33 | 16.54
Summer | -10.22 | 0.74 | 0.42 | -0.31 | 12.19
Winter | -1.71 | 4.69 | 3.80 | -0.45 | 35.87

Finally, estimated parameters of the latent variable and associated
uncertainty are $\hat{\theta}=1.13^{\circ}$ $(std=0.02^{\circ})$ and
$\hat{g}=+0.74$ $(std<0.01)$, which identify the scale parameter of the
exponential covariance function and the temporal parameter of the Markovian
process, respectively. Note that all stations share the same $g$ coefficient
and the latent process is stationary. The magnitude of $g$ indicates that the
time series has components at low frequency, meaning that the changes over
time are relatively smooth.

From residual diagnostics analysis, in Figure 4a the distribution of the
residuals $\epsilon(s,t)$ grouped by the station are shown through boxplots.
Their distributions appear centred around zero. On the other hand, looking at
the residual $\epsilon(s,t)$ globally, shown in Figure 4b, the distribution is
non-Gaussian. In particular, kurtosis is larger than three ($ku=10.4$) and
skewness greater than 0 ($sk=0.14$). The residuals are further investigated
through the ACF, and we note that there is a negative autocorrelation at lag-1
of $-0.2$. Finally, we calculate the daily $RMSE_{t}$ using Eq. 4. As
expected, we note that there is a seasonal component in the residuals. In
particular, the variance is higher in the winter than in summer. To obtain the
estimate $\hat{\sigma}_{\epsilon}(t)$ we use an interpolation spline over the
monthly averages of daily $RMSE_{t}$. Considering that autocorrelation is
negative, the $MSE$ upper bound of Eq. 6 holds true. More details on ACF
computed on the residuals are reported in the Appendix C.

(a)

(b)

Figure 4: Residual $\epsilon(\bm{s},t)$ of Eq. 1 grouped by station (Figure
4a) and globally (Figure 4b). The outliers (red cross) are identified by the
Whisker length based on the first and third quartiles.  Figure 5: Daily
$RMSE_{t}$ (dotted blue line) computed by Eq. 4. The daily estimate is then
filtered by MA(p = 30) (dotted red line), monthly average step function (pink
cross) and by an interpolation spline over the monthly averages (dotted purple
line).

###  3.4 Scenario Analysis

We consider the "What-if" scenario analysis approach: "What would happen if
the NH3 reductions were fully implemented at time 0 (1st January 2016)?".
Considering those $NH_{3}$ emissions and the AQ stations are concentrated in
the Po Valley area as shown in Figure 2, we evaluate the scenario analysis
only in areas under the 640 $[m]$ altitude and not classified as a forest.
Table 1 summarises the proposed scenarios. For each scenario, we map the
average and the uncertainty of the PM2.5 changes due to NH3 reductions on a
regular grid of 0.1∘ $\times$ 0.1∘ over the studied area. Maps are built by
Eq. 6. We focus on the winter period (452 daily predictions for each pixel),
which is characterised by high PM2.5 concentration (see: Table 4) and the
biggest effect of NH3 (see Table 5). Note that the average is not spread over
the mountain area where the magnitude of the $NH_{3}$ emissions is low.

Figure 6a depicts the average PM2.5 changes for the PRIA scenario while Figure
6b shows the associate uncertainty. In an analogous way, Figure 7a and Figure
7b show the average and the associate uncertainty of the PM2.5 changes for the
Strong scenario. Both scenarios show that the main reduction effect is
obtained in the southeast area of the region where the $NH_{3}$ emissions are
the highest. In particular, in some areas, the reduction of PM2.5
concentrations reaches more than 7 $\mu g/m^{3}$. For example, in Brescia or
Cremona, located in the southeast part of the region, the average of the PM2.5
concentrations are ones of the highest of the region, and an average reduction
of 6-7 $\mu g/m^{3}$ in winter period can significantly improve the local air
quality, reducing the risks associated to the health of the population. On the
other side, heavily urbanised land such as the metropolitan area of Milan, do
not benefit from the reduction of ammonia emissions as they don’t present high
levels of NH3 emissions. A more marked reduction of NH3 emissions (Strong
scenario) leads to more significant improvements in air quality, achieving
local PM2.5 reductions close to 15 $\mu g/m^{3}$.

(a)

(b)

Figure 6: PRIA Scenario (-26%): Panel (6a): $\bar{\Delta}_{Y}$ computed in the
winter season (452 daily observations for each pixel). Panel (6b): uncertainty
associated with the average reductions ($RMSE[\bar{\Delta}_{Y}]$) The pink
stars depict the provincial capitals. The model is only evaluated below the
elevation of 640 $[m]$ in a non-forested area.

(a)

(b)

Figure 7: Strong Scenario (-50%). Panel (7a): $\bar{\Delta}_{Y}$ computed in
the winter season (452 daily observations for each pixel). Panel (7b):
uncertainty associated with the average reductions ($RMSE[\bar{\Delta}_{Y}]$).
The pink stars depict the provincial capitals. The model is only evaluated
below the elevation of 640 $[m]$ in a non-forested area.

Furthermore, we consider the distribution of the daily PM2.5 changes during
winter time aggregated by provinces and land type. Fig. 8 summarises the
boxplot of the daily $\Delta_{Y}$, for the PRIA scenario while Figure 9
summarises the boxplot for the Strong scenario. The main reduction effects are
obtained in Brescia, Cremona and Mantova provinces which are characterised by
large rural areas, which is coherent with the land type aggregation.

(a)

(b)

Figure 8: Scenario PRIA (-26%). Boxplot of daily ${\Delta}_{Y}$ in the winter
season (452 daily observations for each pixel): grouped by Lombardy provinces
8a; grouped by land type 8b. The number over each boxplot shows the number of
pixels of each category. The model is evaluated only in that area under the
640 $[m]$ altitude which is not a forest.

(a)

(b)

Figure 9: Scenario Strong (-50%). Boxplot of daily ${\Delta}_{Y}$ in the
winter season (452 daily observations for each pixel): grouped by Lombardy
provinces 9a; grouped by land type 9b. The number over each boxplot shows the
number of pixels of each category. The model is evaluated only in that area
under the 640 $[m]$ altitude which is not a forest.

Table 7 summarises the averages of PM2.5 reductions aggregated by province and
the Lombardy overall reduction in $\mu g/m^{3}$. In particular, it shows for
each province in the Po Valley basin the number of pixels involved in the
computation ($|D|$), the average prediction of $PM_{2.5}$, $\hat{y}(\bm{s},t)$
and in the last two columns the average PM2.5 changes with the associated
uncertainty in bracket: $\bar{\Delta}_{y}$ and $(MSE[\bar{\Delta}_{y}])$
computed by the Eq. 6. Note that the scenario analysis covers a limited
surface of the Lecco, Varese, and Como provinces. Finally, Table 8 shows the
reduction aggregate by season. As expected, the main effect is obtained in the
winter season, where $NH_{3}$ plays a limited role in the concentration of
$PM_{2.5}$. In the wintertime when the NH3 is less abundant, the overall
estimated reduction over the Po Valley areas for the PRIA scenario is $2.09$
($std=0.10$) while for the Strong scenario is $4.02$ ($std=0.19$). These
results can be compared with the prediction average $39.92$ $\mu g/m^{3}$. So,
considering the PRIA scenario, the PM2.5 reduction is close to 5% while
considering the Strong scenario, the PM2.5 reduction is close to 10%. The
reduction could be achieved through the containment of ammonia emissions
deriving from (i) livestock housing structures, (ii) storage of livestock
manure (iii) distribution of livestock effluents and, (iv) promotion of
conservation agriculture.

Table 7: Winter change aggregate by province in $\mu g/m^{3}$. $|D|$ is the number of pixels involved in the computation, $\hat{y}(\bm{s},t)$ is the average prediction of the $PM_{2.5}$ while the last two columns show change and its uncertainty: $\bar{\Delta}_{y}$ ($RMSE[\bar{\Delta}_{y}]$) for both scenarios. |  | $\hat{y}(\bm{s},t)$ | PRIA (-26%) | Strong (-50%)
---|---|---|---|---
Prov. | $|D|$ | $\mu g/m^{3}$ | $\mu g/m^{3}$ | $\mu g/m^{3}$
VA | 5 | 33.30 | -0.28 (0.08) | -0.54 (0.08)
CO | 4 | 37.46 | -0.37 (0.09) | -0.72 (0.10)
MI | 19 | 37.61 | -0.68 (0.05) | -1.31 (0.07)
BG | 8 | 41.33 | -2.35 (0.13) | -4.52 (0.22)
BS | 21 | 49.32 | -4.56 (0.22) | -8.78 (0.42)
PV | 30 | 35.47 | -0.73 (0.05) | -1.40 (0.07)
CR | 22 | 41.46 | -2.96 (0.15) | -5.70 (0.27)
MN | 24 | 40.55 | -2.87 (0.14) | -5.52 (0.26)
LC | 2 | 34.49 | -0.28 (0.13) | -0.54 (0.13)
LO | 10 | 37.91 | -2.28 (0.12) | -4.39 (0.21)
MB | 4 | 38.34 | -0.35 (0.09) | -0.67 (0.10)
Lombardy | 149 | 39.92 | -2.09 (0.10) | -4.02 (0.19)
Table 8: $\bar{\Delta}_{y}(\bm{s},t)$ aggregate by season and the yearly reduction in $\mu g/m^{3}$. The number of pixels involved in the computation is 149. $\hat{y}(\bm{s},t)$ is the average prediction of the $PM_{2.5}$ while the last two columns show the change and its uncertainty: $\bar{\Delta}_{y}$ ($RMSE[\bar{\Delta}_{y}]$) and the (%) of reductions for both scenarios. | $\hat{y}(\bm{s},t)$ | PRIA(-26%) | Strong(-50%)
---|---|---|---
| $\mu g/m^{3}$ | $\mu g/m^{3}$ | % | $\mu g/m^{3}$ | %
Aut. | 24.54 | -0.61 (0.10) | -2.49% | -1.18 (0.19) | -4.79%
Spr. | 18.53 | 0.02 (0.10) | 0.11% | 0.04 (0.20) | 0.22%
Sum. | 14.00 | -0.33 (0.11) | -2.32% | -0.63 (0.20) | -4.47%
Win. | 39.92 | -2.09 (0.10) | -5.23% | -4.02 (0.19) | -10.06%
Year | 24.18 | -0.75 (0.05) | -3.08% | -1.43 (0.11) | -5.93%

##  4 Conclusion

This study carried out as part of the Agrimonia project, quantified the
reduction in $PM_{2.5}$ concentrations achieved through a reduction in
$NH_{3}$ emissions for the southern area of the Lombardy region. If we look at
the entire area considered and the entire period, against an ammonia emissions
reduction of 26% (resulting from the application of the actions identified in
the regional air quality plan) we obtain a fine particulate reduction of 3%
while a reduction of 50% in ammonia emissions improve the air quality by a 6%.
On the other hand, the absolute reduction in $PM_{2.5}$ depends on the
fraction of fine-particulate mass that is directly ammonia-sensitive. In the
Po basin, the composition of PM2.5 in winter has a very significant secondary
component, usually exceeding 50%. This percentage can even be higher during
accumulation episodes under unfavourable meteorological conditions for
dispersion. As a result, the greatest reductions are obtained during the
winter period (over 5% and 10% for both scenarios, respectively). Looking at
the spatial distribution of the reductions, it can be seen that the most
benefited areas are enclosed in the southeastern area of the region, those
areas corresponding to the high density of livestock farms and therefore of
ammonia emissions. Consequently, Lombardy has the largest potential of
reducing winter averaged PM2.5 by strongly controlling NH3 emissions which are
considered beneficial to human health.

##  Appendix A Derivation of MSE upper bound

Starting from Eq. 5 of the daily changes, if we assume $\epsilon_{t}$ and
$\epsilon_{t}^{\alpha}$ are independent with the same variance, then:

$MSE[\Delta_{y}(\bm{s},t)]=\Delta_{\bm{x}}^{\prime}\bm{\Sigma}_{\hat{\bm{\beta}}}\Delta_{\bm{x}}+2V[\epsilon_{t}]$
(9)

where $\epsilon_{t}\sim N(0,\sigma_{\epsilon}^{2}(t))$.

We are interested in aggregating the daily changes over space subset
$\mathcal{D}$ and time interval $\mathcal{I}$. In this setting, we compute the
variance of the average $\bar{\epsilon}$ as follows:

$V[\bar{\epsilon}]=\frac{1}{|D|^{2}|I|^{2}}\frac{1+\gamma}{1-\gamma}\sum_{\bm{s}\in\mathcal{D}}\sum_{t\in\mathcal{I}}\hat{\sigma}_{\epsilon}^{2}(s,t)$
(10)

where $\hat{\sigma}_{\epsilon}^{2}(t)$ is an estimate of the measurement error
variance; $\gamma$ is the autocorrelation coefficients. If the errors are IID
then the $\gamma=0$ otherwise if there is autocorrelation in the residual we
can adjust the variance by setting $\gamma$ accordingly. In particular, if the
autocorrelation is negative, we obtain the upper bound in Eq. 6 by setting
$\gamma=0$.

##  Appendix B Temporal dynamics PM2.5

The temporal dynamics of the response variable, PM2.5 are investigated through
the Partial Auto-Correlation Function (PACF). The PACF shown in Figure 10
highlights the first lag (1 day) as the most significant autocorrelated
component. We compute the PACF for each station and then we summarise the
partial autocorrelation coefficients for all stations by using the boxplot
representation. For all lags, the confidence bounds are computed by $0\pm
2\tilde{\sigma}$ where $\tilde{\sigma}$ is the estimated standard error of the
sample partial autocorrelation.

Figure 10: Partial Auto-Correlation Function (PACF) calculated by stations for
the first 60 lags (days) and summarised through box plots. The i-th box plot
summarises the partial autocorrelation calculated for each station at i-th
lag.

##  Appendix C Autocorrelation function

To check the autocorrelation in the residuals, we compute the Autocorrelation
Function (ACF). The ACF is shown in Figure 11. Since we have residuals for
each station, we compute the ACF for each station and then we summarise the
autocorrelation coefficients for all stations by using the boxplot
representation.

Figure 11: Autocorrelation Function (ACF) computed by stations residual for
the first 60 lags (days) and summarised through box plots. The i-th box plot
summarises the autocorrelation calculated for each station at i-th lag.

## Credit author statement

Alessandro Fssò: Conceptualisation, Methodology, Supervision. Jacopo
Rodeschini: Formal analysis, Software, Data curation. Alessandro Fusta Moro:
Data curation, Revision of the article. Francesco Finazzi: Revision of the
article.

## Acknowledgement

This research was co-funded by Fondazione Cariplo under the grant 2020–4066
“AgrImOnIA: the impact of agriculture on air quality and the COVID-19
pandemic” from the “Data Science for Science and Society” program and by the
European Union - NextGenerationEU, in the framework of the “GRINS - Growing
Resilient, INclusive and Sustainable” project (GRINS PE00000018 – CUP
F83C22001720001). The views and opinions expressed are solely those of the
authors and do not necessarily reflect those of the European Union, nor can
the European Union be held responsible for them.

## References

  * Beven [2019] K. Beven.  Validation and equifinality.  _Computer Simulation Validation: Fundamental Concepts, Methodological Frameworks, and Philosophical Perspectives_ , pages 791–809, 2019\.
  * Calculli et al. [2015] C. Calculli, A. Fassò, F. Finazzi, A. Pollice, and A. Turnone.  Maximum likelihood estimation of the multivariate hidden dynamic geostatistical model with application to air quality in apulia, italy.  _Environmetrics_ , 26(6):406–417, 2015.
  * Clappier et al. [2021] A. Clappier, P. Thunis, M. Beekmann, J. Putaud, and A. de Meij.  Impact of sox, nox and nh3 emission reductions on pm2. 5 concentrations across europe: Hints for future measure development.  _Environment International_ , 156:106699, 2021.
  * Cressie and Wikle [2015] N. Cressie and C. K. Wikle.  _Statistics for spatio-temporal data_.  John Wiley & Sons, 2015.
  * Cressie et al. [2010] N. Cressie, T. Shi, and E. L. Kang.  Fixed rank filtering for spatio-temporal data.  _Journal of Computational and Graphical Statistics_ , 19(3):724–745, 2010.
  * De Meij et al. [2009] A. De Meij, P. Thunis, B. Bessagnet, and C. Cuvelier.  The sensitivity of the chimere model to emissions reduction scenarios on air quality in northern italy.  _Atmospheric environment_ , 43(11):1897–1907, 2009.
  * Durbin and Koopman [2012] J. Durbin and S. J. Koopman.  _Time Series Analysis by State Space Methods_ , volume 38.  OUP Oxford, 2012.
  * European Environment Agency (2022) [EEA] European Environment Agency (EEA).  Premature deaths due to exposure to fine particulate matter in Europe (8th EAP), 2022.  URL https://www.eea.europa.eu/ims/health-impacts-of-exposure-to.
  * Fassò [2023] A. Fassò.  To what extent airborne particulate matters are influenced by ammonia and nitrogen oxides?  _arXiv preprint arXiv:2310.09302_ , 2023.
  * Fassò et al. [2023a] A. Fassò, J. Rodeschini, A. Fusta Moro, Q. Shaboviq, P. Maranzano, M. Cameletti, F. Finazzi, N. Golini, R. Ignaccolo, and P. Otto.  AgrImOnIA: Open Access dataset correlating livestock and air quality in the Lombardy region, Italy, May 2023a.  URL https://doi.org/10.5281/zenodo.7956006.
  * Fassò et al. [2023b] A. Fassò, J. Rodeschini, A. F. Moro, Q. Shaboviq, P. Maranzano, M. Cameletti, F. Finazzi, N. Golini, R. Ignaccolo, and P. Otto.  Agrimonia: a dataset on livestock, meteorology and air quality in the lombardy region, italy.  _Scientific Data_ , 2023b.  10.1038/s41597-023-02034-0.
  * Inness et al. [2019] A. Inness, M. Ades, A. Agustí-Panareda, J. Barré, A. Benedictow, A.-M. Blechschmidt, J. J. Dominguez, R. Engelen, H. Eskes, J. Flemming, V. Huijnen, L. Jones, Z. Kipling, S. Massart, M. Parrington, V.-H. Peuch, M. Razinger, S. Remy, M. Schulz, and M. Suttie.  The CAMS reanalysis of atmospheric composition.  _Atmos. Chem. Phys._ , 19(6):3515–3556, 2019\.  https://doi.org/10.5194/acp-19-3515-2019.  URL https://acp.copernicus.org/articles/19/3515/2019/.
  * Lonati and Cernuschi [2020] G. Lonati and S. Cernuschi.  Temporal and spatial variability of atmospheric ammonia in the lombardy region (northern italy).  _Atmospheric Pollution Research_ , 11(12):2154–2163, 2020.
  * Lovarelli et al. [2021] D. Lovarelli, D. Fugazza, M. Costantini, C. Conti, G. Diolaiuti, and M. Guarino.  Comparison of ammonia air concentration before and during the spread of covid-19 in lombardy (italy) using ground-based and satellite data.  _Atmospheric Environment_ , 259:118534, 2021.
  * Mannarini et al. [2022] G. Mannarini, M. L. Salinas, L. Carelli, and A. Fassò.  How covid-19 affected ghg emissions of ferries in europe.  _Sustainability_ , 14(9):5287, 2022.
  * Megaritis et al. [2013] A. Megaritis, C. Fountoukis, P. Charalampidis, C. Pilinis, and S. N. Pandis.  Response of fine particulate matter concentrations to changes of emissions and temperature in europe.  _Atmospheric Chemistry and Physics_ , 13(6):3423–3443, 2013.
  * Nowak and Welsh [2020] G. Nowak and A. Welsh.  Improved prediction for a spatio-temporal model.  _Environmental and Ecological Statistics_ , 27:631–648, 2020.
  * Otto et al. [2023] P. Otto, A. F. Moro, J. Rodeschini, Q. Shaboviq, R. Ignaccolo, N. Golini, M. Cameletti, P. Maranzano, F. Finazzi, and A. Fassò.  Spatiotemporal modelling of pm2.5 concentrations in lombardy (italy) – a comparative study, 2023.
  * Pozzer et al. [2017] A. Pozzer, A. P. Tsimpidi, V. A. Karydis, A. De Meij, and J. Lelieveld.  Impact of agricultural emission reductions on fine-particulate matter and public health.  _Atmospheric Chemistry and Physics_ , 17(20):12813–12826, 2017.
  * Rue et al. [2009] H. Rue, S. Martino, and N. Chopin.  Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 71(2):319–392, 2009.
  * Veratti et al. [2023] G. Veratti, M. Stortini, R. Amorati, L. Bressan, G. Giovannini, S. Bande, F. Bissardella, S. Ghigo, E. Angelino, L. Colombo, et al.  Impact of nox and nh3 emission reduction on particulate matter across po valley: A life-ip-prepair study.  _Atmosphere_ , 14(5):762, 2023.
  * Verweij et al. [2023] R. W. Verweij, E. Pisoni, A. van Pul, P. Thunis, and E. van der Swaluw.  A fast method for ammonia emission reduction scenario studies: Combining emep4nl with the sherpa tool.  _Atmospheric Environment_ , 306:119805, 2023.
  * Wang et al. [2021] Y. Wang, F. Finazzi, and A. Fassò.  D-stem v2: A software for modeling functional spatio-temporal data.  _Journal of Statistical Software_ , 99(10):1–29, 2021.  10.18637/jss.v099.i10.  URL https://www.jstatsoft.org/index.php/jss/article/view/v099i10.
  * Wikle [2003] C. K. Wikle.  Hierarchical models in environmental science.  _International Statistical Review_ , 71(2):181–199, 2003.


# Shallow Syntax in Deep Water

Swabha Swayamdipta♠ Matthew Peters♣
Brendan Roof♣ Chris Dyer♡ Noah A. Smith♢♣

♠Language Technologies Institute, Carnegie Mellon University
♣Allen Institute for Artificial Intelligence
♢Paul G. Allen School of Computer Science & Engineering, University of
Washington
♡Google DeepMind
<EMAIL_ADDRESS>Work done during an internship at the
Allen Institute for Artificial Intelligence.

###### Abstract

Shallow syntax provides an approximation of phrase-syntactic structure of
sentences; it can be produced with high accuracy, and is computationally cheap
to obtain. We investigate the role of shallow syntax-aware representations for
NLP tasks using two techniques. First, we enhance the ELMo architecture
(Peters et al., 2018b) to allow pretraining on predicted shallow syntactic
parses, instead of just raw text, so that contextual embeddings make use of
shallow syntactic context. Our second method involves shallow syntactic
features obtained automatically on downstream task data. Neither approach
leads to a significant gain on any of the four downstream tasks we considered
relative to ELMo-only baselines. Further analysis using black-box probes from
Liu et al. (2019) confirms that our shallow-syntax-aware contextual embeddings
do not transfer to linguistic tasks any more easily than ELMo’s embeddings. We
take these findings as evidence that ELMo-style pretraining discovers
representations which make additional awareness of shallow syntax redundant.

##  1 Introduction

The NLP community is revisiting the role of linguistic structure in
applications with the advent of contextual word representations (cwrs) derived
from pretraining language models on large corpora Peters et al. (2018a);
Radford et al. (2018); Howard and Ruder (2018); Devlin et al. (2018). Recent
work has shown that downstream task performance may benefit from explicitly
injecting a syntactic inductive bias into model architectures (Kuncoro et al.,
2018), even when cwrs are also used Strubell et al. (2018). However, high
quality linguistic structure annotation at a large scale remains expensive—a
trade-off needs to be made between the quality of the annotations and the
computational expense of obtaining them. Shallow syntactic structures (Abney,
1991; also called chunk sequences) offer a viable middle ground, by providing
a flat, non-hierarchical approximation to phrase-syntactic trees (see Fig. 1
for an example). These structures can be obtained efficiently, and with high
accuracy, using sequence labelers. In this paper we consider shallow syntax to
be a proxy for linguistic structure.

Figure 1:  A sentence with its phrase-syntactic tree (brown) and shallow
syntactic (chunk) annotations (red). Nodes in the tree which percolate down as
chunk labels are in red. Not all tokens in the sentence get chunk labels;
e.g., punctuation is not part of a chunk.

While shallow syntactic chunks are almost as ubiquitous as part-of-speech tags
in standard NLP pipelines Jurafsky and Martin (2000), their relative merits in
the presence of cwrs remain unclear. We investigate the role of these
structures using two methods. First, we enhance the ELMo architecture (Peters
et al., 2018b) to allow pretraining on predicted shallow syntactic parses,
instead of just raw text, so that contextual embeddings make use of shallow
syntactic context (§2). Our second method involves classical addition of chunk
features to cwr-infused architectures for four different downstream tasks
(§3). Shallow syntactic information is obtained automatically using a highly
accurate model (97% $F_{1}$ on standard benchmarks). In both settings, we
observe only modest gains on three of the four downstream tasks relative to
ELMo-only baselines (§4).

Recent work has probed the knowledge encoded in cwrs and found they capture a
surprisingly large amount of syntax Blevins et al. (2018); Liu et al. (2019);
Tenney et al. (2019). We further examine the contextual embeddings obtained
from the enhanced architecture and a shallow syntactic context, using black-
box probes from Liu et al. (2019). Our analysis indicates that our shallow-
syntax-aware contextual embeddings do not transfer to linguistic tasks any
more easily than ELMo embeddings (§4.2).

Overall, our findings show that while shallow syntax can be somewhat useful,
ELMo-style pretraining discovers representations which make additional
awareness of shallow syntax largely redundant.

##  2 Pretraining with Shallow Syntactic Annotations

We briefly review the shallow syntactic structures used in this work, and then
present a model architecture to obtain embeddings from shallow Syntactic
Context (mSynC).

###  2.1 Shallow Syntax

Base phrase chunking is a cheap sequence-labeling–based alternative to full
syntactic parsing, where the sequence consists of non-overlapping labeled
segments (Fig. 1 includes an example.) Full syntactic trees can be converted
into such shallow syntactic chunk sequences using a deterministic procedure
Jurafsky and Martin (2000). Tjong Kim Sang and Buchholz (2000) offered a rule-
based transformation deriving non-overlapping chunks from phrase-structure
trees as found in the Penn Treebank (Marcus et al., 1993). The procedure
percolates some syntactic phrase nodes from a phrase-syntactic tree to the
phrase in the leaves of the tree. All overlapping embedded phrases are then
removed, and the remainder of the phrase gets the percolated label—this
usually corresponds to the head word of the phrase.

In order to obtain shallow syntactic annotations on a large corpus, we train a
BiLSTM-CRF model Lample et al. (2016); Peters et al. (2017), which achieves
97% $F_{1}$ on the CoNLL 2000 benchmark test set. The training data is
obtained from the CoNLL 2000 shared task Tjong Kim Sang and Buchholz (2000),
as well as the remaining sections (except §23 and §20) of the Penn Treebank,
using the official script for chunk
generation.111https://www.clips.uantwerpen.be/conll2000/chunking/ The standard
task definition from the shared task includes eleven chunk labels, as shown in
Table 1.

Label | % Occurrence
---|---
Noun Phrase (NP) | 51.7
Verb Phrase (VP) | 20.0
Prepositional Phrase (PP) | 19.8
Adverbial Phrase (ADVP) | 3.7
Subordinate Clause (SBAR) | 2.1
Adjective Phrase (ADJP) | 1.9
Verb Particles (PRT) | 0.5
Conjunctive Phrase (CONJ) | 0.06
Interjective Phrase (INTJ) | 0.03
List Marker (LST) | 0.01
Unlike Coordination Phrase (UCP) | 0.002
Table 1: Shallow syntactic chunk phrase types from CoNLL 2000 shared task
(Tjong Kim Sang and Buchholz, 2000) and their occurrence % in the training
data.

###  2.2 Pretraining Objective

Traditional language models are estimated to maximize the likelihood of each
word $x_{i}$ given the words that precede it, $p(x_{i}\mid\bm{x}_{<i})$. Given
a corpus that is annotated with shallow syntax, we propose to condition on
both the preceding words _and_ their annotations.

We associate with each word $x_{i}$ three additional variables (denoted
$c_{i}$): the indices of the beginning and end of the last completed chunk
_before_ $x_{i}$, and its label. For example, in Fig. 2, $c_{4}=\langle
3,3,\text{VP}\rangle$ for $x_{4}=\text{the}$. Chunks, $\bm{c}$ are only used
as conditioning context via $p(x_{i}\mid\bm{x}_{<i},\bm{c}_{\leqslant i})$;
they are not predicted.222A different objective could consider predicting the
next chunks, along with the next word. However, this chunker would have access
to strictly less information than usual, since the entire sentence would no
longer be available. Because the $\bm{c}$ labels depend on the entire sentence
through the CRF chunker, conditioning each word’s probability on any
$\bm{c}_{i}$ means that our model is, strictly speaking, not a language model,
and it can no longer be meaningfully evaluated using perplexity.

A right-to-left model is constructed analogously, conditioning on
$\bm{c}_{\geqslant i}$ alongside $\bm{x}_{>i}$. Following Peters et al.
(2018a), we use a joint objective maximizing data likelihood objectives in
both directions, with shared softmax parameters.

###  2.3 Pretraining Model Architecture

Our model uses two encoders: $e_{\mathit{seq}}$ for encoding the sequential
history ($\bm{x}_{<i}$), and $e_{\mathit{syn}}$ for shallow syntactic (chunk)
history ($\bm{c}_{\leqslant i}$). For both, we use transformers Vaswani et al.
(2017), which consist of large feedforward networks equipped with multiheaded
self-attention mechanisms.

Figure 2:  Model architecture for pretraining with shallow syntax. A
sequential encoder converts the raw text into cwrs (shown in blue). Observed
shallow syntactic structure (chunk boundaries and labels, shown in red) are
combined with these cwrs in a shallow syntactic encoder to get contextualized
representations for chunks (shown in orange). Both representations are passed
through a projection layer to get mSynC embeddings (details shown only in some
positions, for clarity), used both for computing the data likelihood, as
shown, as well as in downstream tasks.

As inputs to $e_{\mathit{seq}}$, we use a context-independent embedding,
obtained from a CNN character encoder Kim et al. (2016) for each token
$x_{i}$. The outputs $\bm{h}_{i}$ from $e_{\mathit{seq}}$ represent words in
context.

Next, we build representations for (observed) chunks in the sentence by
concatenating a learned embedding for the chunk label with $\bm{h}$s for the
boundaries and applying a linear projection ($f_{\mathit{proj}}$). The output
from $f_{\mathit{proj}}$ is input to $e_{\mathit{syn}}$, the shallow syntactic
encoder, and results in contextualized chunk representations, $\bm{g}$. Note
that the number of chunks in the sentence is less than or equal to the number
of tokens.

Each $\bm{h}_{i}$ is now concatentated with $\bm{g}_{c_{i}}$, where
$\bm{g}_{c_{i}}$ corresponds to $c_{i}$, the last chunk before position $i$.
Finally, the output is given by
$\mbox{{mSynC}}_{i}={u}_{\mathit{proj}}(\bm{h}_{i},\bm{g}_{c_{i}})=\bm{W}^{\top}[\bm{h}_{i};\bm{g}_{c_{i}}]$,
where $\bm{W}$ is a model parameter. For training, $\mbox{{mSynC}}_{i}$ is
used to compute the probability of the next word, using a sampled softmax
Bengio et al. (2003). For downstream tasks, we use a learned linear weighting
of all layers in the encoders to obtain a task-specific mSynC, following
Peters et al. (2018a).

##### Staged parameter updates

Jointly training both the sequential encoder $e_{\mathit{seq}}$, and the
syntactic encoder $e_{\mathit{syn}}$ can be expensive, due to the large number
of parameters involved. To reduce cost, we initialize our sequential cwrs
$\bm{h}$, using pretrained embeddings from ELMo-transformer. Once initialized
as such, the encoder is fine-tuned to the data likelihood objective (§2.2).
This results in a staged parameter update, which reduces training duration by
a factor of 10 in our experiments. We discuss the empirical effect of this
approach in §4.3.

##  3 Shallow Syntactic Features

Our second approach incorporates shallow syntactic information in downstream
tasks via token-level chunk label embeddings. Task training (and test) data is
automatically chunked, and chunk boundary information is passed into the task
model via BIOUL encoding of the labels. We add randomly initialized chunk
label embeddings to task-specific input encoders, which are then fine-tuned
for task-specific objectives. This approach does not require a shallow
syntactic encoder or chunk annotations for pretraining cwrs, only a chunker.
Hence, this can more directly measure the impact of shallow syntax for a given
task.333In contrast, in §2, the shallow-syntactic encoder itself, as well as
predicted chunk quality on the large pretraining corpus could affect
downstream performance.

##  4 Experiments

Our experiments evaluate the effect of shallow syntax, via contextualization
(mSynC, §2) and features (§3). We provide comparisons with four
baselines—ELMo-transformer Peters et al. (2018b), our reimplementation of the
same, as well as two cwr-free baselines, with and without shallow syntactic
features. Both ELMo-transformer and mSynC are trained on the 1B word benchmark
corpus (Chelba et al., 2013); the latter also employs chunk annotations
(§2.1). Experimental settings are detailed in Appendix §A.1.

###  4.1 Downstream Task Transfer

We employ four tasks to test the impact of shallow syntax. The first three,
namely, coarse and fine-grained named entity recognition (NER), and
constituency parsing, are _span-based_ ; the fourth is a sentence-level
sentiment classification task. Following Peters et al. (2018a), we do not
apply finetuning to task-specific architectures, allowing us to do a
controlled comparison with ELMo. Given an identical base architecture across
models for each task, we can attribute any difference in performance to the
incorporation of shallow syntax or contextualization. Details of downstream
architectures are provided below, and overall dataset statistics for all tasks
is shown in the Appendix, Table 5.

|  | NER
---
| Fine-grained NER
---
| Constituency Parsing
---
| Sentiment
---
Baseline (no cwr) | 88.1 $\pm$ 0.27 | 78.5 $\pm$ 0.19 | 88.9 $\pm$ 0.05 | 51.6 $\pm$ 1.63
\+ shallow syn. features | 88.6 $\pm$ 0.22 | 78.9 $\pm$ 0.13 | 90.8 $\pm$ 0.14 | 51.1 $\pm$ 1.39
ELMo-transformer (Peters et al., 2018b) | 91.1 $\pm$ 0.26 | — | 93.7 $\pm$ 0.00 | —
ELMo-transformer (our reimplementation) | 91.5 $\pm$ 0.25 | 85.7 $\pm$ 0.08 | 94.1 $\pm$ 0.06 | 53.0 $\pm$ 0.72
\+ shallow syn. features | 91.6 $\pm$ 0.40 | 85.9 $\pm$ 0.28 | 94.3 $\pm$ 0.03 | 52.6 $\pm$ 0.54
Shallow syn. contextualization (mSynC) | 91.5 $\pm$ 0.19 | 85.9 $\pm$ 0.20 | 94.1 $\pm$ 0.07 | 53.0 $\pm$ 1.07
Table 2:  Test-set performance of ELMo-transformer Peters et al. (2018b), our reimplementation, and mSynC, compared to baselines without cwr. Evaluation metric is $F_{1}$ for all tasks except sentiment, which reports accuracy. Reported results show the mean and standard deviation across 5 runs for coarse-grained NER and sentiment classification and 3 runs for other tasks.  | CCG |  | PTB
---
POS
| EWT
---
POS
Chunk | NER |  | Sem.
---
Tagging
| Gramm.
---
Err. D
| Prep.
---
Role
| Prep.
---
Func.
| Event
---
Fact.
ELMo-transformer | 92.68 | 97.09 | 95.13 | 92.18 | 81.21 | 93.78 | 30.80 | 72.81 | 82.24 | 70.88
mSynC | 92.03 | 96.91 | 94.64 | 96.89 | 79.98 | 93.03 | 30.86 | 70.83 | 82.67 | 70.39
Table 3: Test performance of ELMo-transformer Peters et al. (2018b) vs. mSynC
on several linguistic probes from Liu et al. (2019). In each case, performance
of the best layer from the architecture is reported. Details on the probes can
be found in §4.2.1.

##### NER

We use the English portion of the CoNLL 2003 dataset (Tjong Kim Sang and De
Meulder, 2003), which provides named entity annotations on newswire data
across four different entity types (PER, LOC, ORG, MISC). A bidirectional
LSTM-CRF architecture (Lample et al., 2016) and a BIOUL tagging scheme were
used.

##### Fine-grained NER

The same architecture and tagging scheme from above is also used to predict
fine-grained entity annotations from OntoNotes 5.0 (Weischedel et al., 2011).
There are 18 fine-grained NER labels in the dataset, including regular named
entitities as well as entities such as date, time and common numerical
entries.

##### Phrase-structure parsing

We use the standard Penn Treebank splits, and adopt the span-based model from
Stern et al. (2017). Following their approach, we used predicted part-of-
speech tags from the Stanford tagger (Toutanova et al., 2003) for training and
testing. About 51% of phrase-syntactic constituents align exactly with the
predicted chunks used, with a majority being single-width noun phrases. Given
that the rule-based procedure used to obtain chunks only propagates the phrase
type to the head-word and removes all overlapping phrases to the right, this
is expected. We did not employ jack-knifing to obtain predicted chunks on PTB
data; as a result there might be differences in the quality of shallow syntax
annotations between the train and test portions of the data.

##### Sentiment analysis

We consider fine-grained (5-class) classification on Stanford Sentiment
Treebank (Socher et al., 2013). The labels are negative, somewhat_negative,
neutral, positive and somewhat_positive. Our model was based on the
biattentive classification network (McCann et al., 2017). We used all phrase
lengths in the dataset for training, but test results are reported only on
full sentences, following prior work.

Results are shown in Table 2. Consistent with previous findings, cwrs offer
large improvements across all tasks. Though helpful to span-level task models
without cwrs, shallow syntactic features offer little to no benefit to ELMo
models. mSynC’s performance is similar. This holds even for phrase-structure
parsing, where (gold) chunks align with syntactic phrases, indicating that
task-relevant signal learned from exposure to shallow syntax is _already_
learned by ELMo. On sentiment classification, chunk features are slightly
harmful on average (but variance is high); mSynC again performs similarly to
ELMo-transformer. Overall, the performance differences across all tasks are
small enough to infer that shallow syntax is not particularly helpful when
using cwrs.

###  4.2 Linguistic Probes

We further analyze whether awareness of shallow syntax carries over to other
linguistic tasks, via probes from Liu et al. (2019). Probes are linear models
trained on frozen cwrs to make predictions about linguistic (syntactic and
semantic) properties of words and phrases. Unlike §4.1, there is minimal
downstream task architecture, bringing into focus the transferability of cwrs,
as opposed to task-specific adaptation.

####  4.2.1 Probing Tasks

The ten different probing tasks we used include CCG supertagging Hockenmaier
and Steedman (2007), part-of-speech tagging from PTB Marcus et al. (1993) and
EWT (Universal Depedencies Silveira et al., 2014), named entity recognition
Tjong Kim Sang and De Meulder (2003), base-phrase chunking Tjong Kim Sang and
Buchholz (2000), grammar error detection Yannakoudakis et al. (2011), semantic
tagging Bjerva et al. (2016), preposition supersense identification Schneider
et al. (2018), and event factuality detection Rudinger et al. (2018). Metrics
and references for each are summarized in Table 6. For more details, please
see Liu et al. (2019).

Results in Table 3 show ten probes. Again, we see the performance of baseline
ELMo-transformer and mSynC are similar, with mSynC doing slightly worse on 7
out of 9 tasks. As we would expect, on the probe for predicting chunk tags,
mSynC achieves 96.9 $F_{1}$ vs. 92.2 $F_{1}$ for ELMo-transformer, indicating
that mSynC is indeed encoding shallow syntax. Overall, the results further
confirm that explicit shallow syntax does not offer any benefits over ELMo-
transformer.

###  4.3 Effect of Training Scheme

We test whether our staged parameter training (§2.3) is a viable alternative
to an end-to-end training of both $e_{\mathit{syn}}$ and $e_{\mathit{seq}}$.
We make a further distinction between fine-tuning $e_{\mathit{seq}}$ vs. not
updating it at all after initialization (frozen).

| Model |  | Fine-grained
---
NER $F_{1}$
end-to-end | ELMo | 86.90 $\pm$ 0.11
mSynC end-to-end | 86.89 $\pm$ 0.04
staged |  mSynC frozen | 87.36 $\pm$ 0.02
mSynC fine-tuned | 87.44 $\pm$ 0.07
Table 4: Validation $F_{1}$ for fine-grained NER across syntactic pretraining
schemes, with mean and standard deviations across 3 runs.

Downstream validation-set $F_{1}$ on fine-grained NER, reported in Table 4,
shows that the end-to-end strategy lags behind the others, perhaps indicating
the need to train longer than 10 epochs. However, a single epoch on the
1B-word benchmark takes 36 hours on 2 Tesla V100s, making this prohibitive.
Interestingly, the frozen strategy, which takes the least amount of time to
converge (24 hours on 1 Tesla V100), also performs almost as well as fine-
tuning.

##  5 Conclusion

We find that exposing cwr-based models to shallow syntax, either through new
cwr learning architectures or explicit pipelined features, has little effect
on their performance, across several tasks. Linguistic probing also shows that
cwrs aware of such structures do not improve task transferability. Our
architecture and methods are general enough to be adapted for richer inductive
biases, such as those given by full syntactic trees (RNNGs; Dyer et al.,
2016), or to different pretraining objectives, such as masked language
modeling (BERT; Devlin et al., 2018); we leave this pursuit to future work.

## References

  * Abney (1991) Steven P Abney. 1991.  Parsing by chunks.  In _Principle-based parsing_ , pages 257–278. Springer.
  * Bengio et al. (2003) Yoshua Bengio, Réjean Ducharme, Pascal Vincent, and Christian Janvin. 2003.  A neural probabilistic language model.  _Journal of Machine Learning Research_ , pages 1137–1155.
  * Bjerva et al. (2016) Johannes Bjerva, Barbara Plank, and Johan Bos. 2016.  Semantic tagging with deep residual networks.  In _Proc. of COLING_.
  * Blevins et al. (2018) Terra Blevins, Omer Levy, and Luke Zettlemoyer. 2018.  Deep rnns encode soft hierarchical syntax.  In _Proc. of ACL_.
  * Chelba et al. (2013) Ciprian Chelba, Tomas Mikolov, Mike Schuster, Qi Ge, Thorsten Brants, Phillipp Koehn, and Tony Robinson. 2013.  One billion word benchmark for measuring progress in statistical language modeling.  ArXiv:1312.3005.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  BERT: pre-training of deep bidirectional transformers for language understanding.  ArXiv:1810.04805.
  * Dyer et al. (2016) Chris Dyer, Adhiguna Kuncoro, Miguel Ballesteros, and Noah A Smith. 2016.  Recurrent neural network grammars.  In _Proc. of NAACL-HLT_.
  * Gardner et al. (2017) Matt Gardner, Joel Grus, Mark Neumann, Oyvind Tafjord, Pradeep Dasigi, Nelson F. Liu, Matthew Peters, Michael Schmitz, and Luke S. Zettlemoyer. 2017\.  AllenNLP: A deep semantic natural language processing platform.  ArXiv:1803.07640.
  * Hockenmaier and Steedman (2007) Julia Hockenmaier and Mark Steedman. 2007.  CCGbank: A corpus of ccg derivations and dependency structures extracted from the Penn Treebank.  _Computational Linguistics_ , 33(3).
  * Howard and Ruder (2018) Jeremy Howard and Sebastian Ruder. 2018.  Universal language model fine-tuning for text classification.  In _Proc. of ACL_.
  * Jurafsky and Martin (2000) Daniel Jurafsky and James H. Martin. 2000.  _Speech and Language Processing: An Introduction to Natural Language Processing, Computational Linguistics, and Speech Recognition_, 1st edition.  Prentice Hall PTR, Upper Saddle River, NJ, USA.
  * Kim et al. (2016) Yoon Kim, Yacine Jernite, David Sontag, and Alexander M. Rush. 2016.  Character-aware neural language models.  In _Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence_ , AAAI’16, pages 2741–2749. AAAI Press.
  * Kuncoro et al. (2018) Adhiguna Kuncoro, Chris Dyer, John Hale, Dani Yogatama, Stephen Clark, and Phil Blunsom. 2018.  Lstms can learn syntax-sensitive dependencies well, but modeling structure makes them better.  In _Proc. of ACL_.
  * Lample et al. (2016) Guillaume Lample, Miguel Ballesteros, Sandeep Subramanian, Kazuya Kawakami, and Chris Dyer. 2016.  Neural architectures for named entity recognition.  In _Proc. of NAACL-HLT_.
  * Liu et al. (2019) Nelson F. Liu, Matt Gardner, Yonatan Belinkov, Matthew Peters, and Noah A. Smith. 2019.  Linguistic knowledge and transferability of contextual representations.  In _Proc. of NAACL-HLT_.
  * Marcus et al. (1993) Mitchell P Marcus, Mary Ann Marcinkiewicz, and Beatrice Santorini. 1993.  Building a large annotated corpus of English: The Penn treebank.  _Computational Linguistics_ , 19(2):313–330.
  * McCann et al. (2017) Bryan McCann, James Bradbury, Caiming Xiong, and Richard Socher. 2017.  Learned in translation: Contextualized word vectors.  In _Proc. of NeurIPS_ , pages 6294–6305.
  * Peters et al. (2017) Matthew Peters, Waleed Ammar, Chandra Bhagavatula, and Russell Power. 2017.  Semi-supervised sequence tagging with bidirectional language models.  In _Proc. of ACL_.
  * Peters et al. (2018a) Matthew Peters, Mark Neumann, Mohit Iyyer, Matt Gardner, Christopher Clark, Kenton Lee, and Luke Zettlemoyer. 2018a.  Deep contextualized word representations.  In _Proc. of NAACL-HLT_.
  * Peters et al. (2018b) Matthew Peters, Mark Neumann, Luke Zettlemoyer, and Wen-tau Yih. 2018b.  Dissecting contextual word embeddings: Architecture and representation.  In _Proc. of EMNLP_ , pages 1499–1509.
  * Radford et al. (2018) Alec Radford, Karthik Narasimhan, Tim Salimans, and Ilya Sutskever. 2018.  Improving language understanding by generative pre-training.
  * Rudinger et al. (2018) Rachel Rudinger, Aaron Steven White, and Benjamin Van Durme. 2018.  Neural models of factuality.  In _Proc. of ACL_.
  * Schneider et al. (2018) Nathan Schneider, Jena D Hwang, Vivek Srikumar, Jakob Prange, Austin Blodgett, Sarah R Moeller, Aviram Stern, Adi Bitan, and Omri Abend. 2018.  Comprehensive supersense disambiguation of English prepositions and possessives.
  * Silveira et al. (2014) Natalia Silveira, Timothy Dozat, Marie-Catherine De Marneffe, Samuel R Bowman, Miriam Connor, John Bauer, and Christopher D Manning. 2014.  A gold standard dependency corpus for English. In _Proc. of LREC_ , pages 2897–2904.
  * Socher et al. (2013) Richard Socher, Alex Perelygin, Jean Wu, Jason Chuang, Christopher D Manning, Andrew Ng, and Christopher Potts. 2013.  Recursive deep models for semantic compositionality over a sentiment treebank.  In _Proc. of EMNLP_.
  * Stern et al. (2017) Mitchell Stern, Jacob Andreas, and Dan Klein. 2017.  A minimal span-based neural constituency parser.  In _Proc. of ACL_.
  * Strubell et al. (2018) Emma Strubell, Patrick Verga, Daniel Andor, David Weiss, and Andrew McCallum. 2018\.  Linguistically-informed self-attention for semantic role labeling.  In _Proc. of EMNLP_.
  * Tenney et al. (2019) Ian Tenney, Patrick Xia, Berlin Chen, Alex Wang, Adam Poliak, R Thomas McCoy, Najoung Kim, Benjamin Van Durme, Sam Bowman, Dipanjan Das, and Ellie Pavlick. 2019\.  What do you learn from context? Probing for sentence structure in contextualized word representations.  In _Proc. of ICLR_.
  * Tjong Kim Sang and Buchholz (2000) Erik F. Tjong Kim Sang and Sabine Buchholz. 2000.  Introduction to the CoNLL-2000 shared task: Chunking.  In _Proc. of CoNLL_.
  * Tjong Kim Sang and De Meulder (2003) Erik F Tjong Kim Sang and Fien De Meulder. 2003.  Introduction to the CoNLL-2003 shared task: Language-independent named entity recognition.  In _Proc. of NAACL_. Association for Computational Linguistics.
  * Toutanova et al. (2003) Kristina Toutanova, Dan Klein, Christopher D. Manning, and Yoram Singer. 2003.  Feature-rich part-of-speech tagging with a cyclic dependency network.  In _Proc. of NAACL_.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Proc. of NeurIPS_ , pages 5998–6008.
  * Weischedel et al. (2013) Ralph Weischedel, Martha Palmer, Mitchell Marcus, Eduard Hovy, Sameer Pradhan, Lance Ramshaw, Nianwen Xue, Ann Taylor, Jeff Kaufman, Michelle Franchini, et al. 2013.  OntoNotes release 5.0 ldc2013t19.  _Linguistic Data Consortium, Philadelphia, PA_.
  * Weischedel et al. (2011) Ralph Weischedel, Sameer Pradhan, Lance Ramshaw, Martha Palmer, Nianwen Xue, Mitchell Marcus, Ann Taylor, Craig Greenberg, Eduard Hovy, Robert Belvin, et al. 2011.  OntoNotes Release 4.0.  _LDC2011T03, Philadelphia, Penn.: Linguistic Data Consortium_.
  * Yannakoudakis et al. (2011) Helen Yannakoudakis, Ted Briscoe, and Ben Medlock. 2011.  A new dataset and method for automatically grading esol texts.  In _Proc. of ACL_.

##  Appendix A Supplemental Material

###  A.1 Hyperparameters

##### ELMo-transformer

Our baseline pretraining model was a reimplementation of that given in Peters
et al. (2018b). Hyperparameters were generally identical, but we trained on
only 2 GPUs with (up to) 4,000 tokens per batch. This difference in batch size
meant we used 6,000 warm up steps with the learning rate schedule of Vaswani
et al. (2017).

#####  mSynC

The function $f_{seq}$ is identical to the 6-layer biLM used in ELMo-
transformer. $f_{syn}$, on the other hand, uses only 2 layers. The learned
embeddings for the chunk labels have 128 dimensions and are concatenated with
the two boundary $\bm{h}$ of dimension 512. Thus $f_{proj}$ maps $1024+128$
dimensions to 512. Further, we did not perform weight averaging over several
checkpoints.

##### Shallow Syntax

The size of the shallow syntactic feature embedding was 50 across all
experiments, initialized uniform randomly.

All model implementations are based on the AllenNLP library Gardner et al.
(2017).

Task | Train | Heldout | Test
---|---|---|---
CoNLL 2003 NER (Tjong Kim Sang and De Meulder, 2003) | 23,499 | 5,942 | 5,648
OntoNotes NER (Weischedel et al., 2013) | 81,828 | 11,066 | 11,257
Penn TreeBank (Marcus et al., 1993) | 39,832 | 1,700 | 2,416
Stanford Sentiment Treebank (Socher et al., 2013) | 8,544 | 1,101 | 2,210
Table 5: Downstream dataset statistics describing the number of train, heldout and test set instances for each task.  Task | Dataset | Metric
---|---|---
CCG Supertagging | CCGBank Hockenmaier and Steedman (2007) | Accuracy
PTB part-of-speech tagging | PennTreeBank Marcus et al. (1993) | Accuracy
EWT part-of-speech tagging | Universal Dependencies Silveira et al. (2014) | Accuracy
Chunking | CoNLL 2000 Tjong Kim Sang and Buchholz (2000) | $F_{1}$
Named Entity Recognition | CoNLL 2003 Tjong Kim Sang and De Meulder (2003) | $F_{1}$
Semantic Tagging | Bjerva et al. (2016) | Accuracy
Grammar Error Detection | First Certificate in English Yannakoudakis et al. (2011) | $F_{1}$
Preposition Supersense Role | STREUSLE 4.0 Schneider et al. (2018) | Accuracy
Preposition Supersense Function | STREUSLE 4.0 Schneider et al. (2018) | Accuracy
Event Factuality Detection | UDS It Happened v2 Rudinger et al. (2018) | Pearson R
Table 6: Dataset and metrics for each probing task from Liu et al. (2019),
corresponding to Table 3.