Deep Hurdle Networks for Zero-Inflated Multi-Target Regression: Application to Multiple Species Abundance Estimation

1 Introduction

Since the Industrial Revolution there has been an increase in biodiversity loss, due to a combination of factors such as agriculture, urbanization, and deforestation, as well as climate change and human introduction of non-native species to ecosystems. Biodiversity loss is a great challenge for humanity, given the importance of biodiversity for sustaining ecosystem services. For example, bird species play a key role in regulating ecosystems by controlling pests, pollinating flowers, spreading seeds and regenerating forests. In a current study it was shown that the bird population in the United States and Canada has fallen by an estimated 29% since 1970 [Rosenberg et al., 2019]. The biomass of top marine predators has declined substantially in many cases, [McCauley et al., 2015], and many marine species are shifting rapidly to new regions in response to changing ocean conditions [Pinsky et al., 2020]. More generally, a recent report from the United Nations warned that about a million animal and plant species face extinction [Brondizio et al., 2019].

∗Contact Author

To protect and conserve species, a key question in computational sustainability concerns understanding the distribution of species across landscapes over time, which gives rise to challenging large-scale spatial and temporal modeling and prediction problems [Gomes et al., 2019]. In particular, ecologists are interested in understanding how species interact with the environment as well as how species interact with each other. Joint species distribution modeling is therefore a computational challenge, as we are interested in simultaneously modeling the correlated distributions of potentially hundreds of species, rather than a single species at a time as traditionally done. Another challenge in joint species distribution modeling is that often outcomes of interest, such as local species abundance in terms of counts or biomass, are sparsely observed, leading to zero-inflated data [Morley et al., 2018]. Zero-inflated data are frequent in many other settings beyond ecology, for example in research studies concerning public health when counting the number of vaccine adverse events [Rose et al., 2006] and prediction of sparse user-item consumption rates [Lichman and Smyth, 2018].

Herein we propose general models for jointly estimating counts or abundance for multiple entities, which are referred to as zero-inflated multi-target regression. While our models are general, we focus on computational sustainability applications, in particular the joint estimation of counts of multiple bird species and biomass of multiple fish species. As discussed above, many of these species have suffered dramatic reductions or changes in geographic distributions in recent years.

Our contributions are:

1. We propose a deep generalization of the conventional model for multi-target regression that simultaneously models zero-inflated data and the correlations among the multiple response variables.
2. We provide an efficient end-to-end learning framework for our model that can be implemented on GPUs.
3. We evaluate our model as well as state-of-the art joint distribution models on two datasets from computational sustainability that concern the distributions of fish and bird species.

2 RelatedWork

A popular model for zero-inflated data is the so-called hurdle model [Mullahy, 1986] in which a Bernoulli distribution governs the binary outcome of whether a response variable has a zero or positive realization; if the realization is positive, the "hurdle" is crossed, and the conditional distribution of the positives is governed by a zero-truncated parametric distribution. Although the hurdle model is popular for zero-inflated data, it has several limitations: 1) the two components of the model are assumed to be independent; 2) it does not explicitly capture the relationships between multiple response variables; and 3) its traditional implementation does not scale. In our work, we address these three key limitations.

A closely related alternative to handle zero-inflated count data is the family of zero-inflated models [Lichman and Smyth, 2018], in which the response variable is also modeled as a mixture of a Bernoulli distribution and a parametric distribution on non-negative integers such as the Poisson or negative binomial distributions. Different from hurdle model, the conditional distribution of a zero-inflated model is not required to be zero-truncated. In other words, while the hurdle model assumes zeros only come from the Bernoulli distribution, a zero-inflated model assumes zeros could come from both the Bernoulli distribution and the conditional distribution. To account for the inherent correlation of response variables, a class of multi-level zero-inflated regression models was presented [Almasi et al., 2016]. With such a model, response variables are organized in a hierarchical structure. Response variables are taken to be independent between clusters. Cluster level and within-cluster correlations of response variables are modeled explicitly through random effects attached to linear predictors. However, characterising a suitable hierarchical structure is nontrivial, and these models also do not explicitly capture covariance relations between response variables.

The general problem of multiple-target regression has been extensively studied [Borchani et al., 2015; Xi et al., 2018; Xu et al., 2019]. Existing methods for multi-output regression can be categorized as: (1) problem transformation methods that transform the multi-output problem into independent single-output problems each solved using a single-output regression algorithm, such as multi-target regressor stacking (MTRS) [Spyromitros-Xioufis et al., 2016] and regressor chains (RC) [Melki et al., 2017]; and (2) algorithm adaptation methods that adapt a specific single-output method to directly handle multi-output dataset, such as multi-objective random forest (MORF) [Kocev et al., 2007], random linear target combination (RLTC) [Tsoumakas et al., 2014], and multi-output support vector regression (MOSVR) [Zhu and Gao, 2018]. An empirical comparison of three representative state-of-the-art multi-output regression learning algorithms, MTRS, RLTC and MORF, is presented in [Xi et al., 2018]. Although these advanced multiple-output regression algorithms exploit some correlation between response variables to improve predictive performance, they do not fully model the covariance structure of the response variables and do not consider zero-inflated data. Our work explicitly models zero-inflated data in multi-target regression and it also explicitly models the covariance underlying the phenomena to better characterize the correlation among entities.

3 Preliminaries

In this section we give a brief introduction to the hurdle model and multivariate probit model. In the following we use R, R0, R+, N, N0 and N+ to denote the reals, nonnegative reals, positive reals, integers, nonnegative integers and positive integers, respectively. If an L-dimensional vector y is integral, we write y ∈ NL. We write y ∈ RL if y is real. Two vectors x, y ∈ RL satisfy x ≼ y (or x ≽ y) iff xj ≤ yj (or xj ≥ yj) for 1 ≤ j ≤ L. Since we consider the problem of species abundance estimation, in the following all label data are assumed to be nonnegative. We denote the probability density function (PDF) and the cumulative density function (CDF) of a multivariate normal distribution N(µ, Σ) as φ(x | µ, Σ) and Φ(x | µ, Σ), respectively.

3.1 Hurdle Model

The hurdle model aims to fit data with two independent distributions: (1) a Bernoulli distribution which governs the binary outcome of a response variable being zero or positive; and (2) a zero-truncated distribution of the positive response variable. Specifically, given a dataset D = {(x, y)(i) | i = 1, . . . , N}, where x ∈ RM is the feature data and y ∈ RL0 (or y ∈ NL0 ) is the label data, and let pj be the probability of yj (1 ≤ j ≤ L) being positive, then we have y′j ∼ Bernoulli(pj), where y′j = 1 if yj > 0 and 0 otherwise. Let f(yj | yj > 0) be the PDF of the zero-truncated distribution of yj. The likelihood of yj is given as L(yj) = Pr(y′j = 1)f(yj | yj > 0).

3.2 Multivariate Probit Model

Given a dataset D = {(x, y′ )(i) | i = 1, . . . , N}, where x ∈ RM is the feature data and y′ ∈ {0, 1}L is the present/absent label data, the multivariate probit model (MVP) [Chen et al., 2018] maps the Bernoulli distribution of each binary response variable y′j to a latent variable rj ∈ R through threshold 0, where r = (r1, . . . , rL) is subject to a multivariate normal distribution:

where r ∼N(µ,Σ). The joint likelihood of y′ is given as

where Aj = (-∞, 0] if y′j = 0, Aj = [0,∞) if y′j = 1. Although there is no closed-form expression for the CDF

of a general multivariate normal distribution, [Chen et al., 2018] proposed an efficient parallel sampling algorithm to approximate it. We can first translate equation (3) into the CDF of a multivariate normal distribution using the affine transformation:

where µ′ = Uµ, Σ′ = UΣU, and U = diag(2y′ - 1) ∈ {-1, 0, 1}L×L is the diagonal matrix with vector 2y′ - 1 as its diagonal. By decomposing Σ′ into I + Σ′′ , where I is the identity matrix, a random variable r ∼ N(0, Σ′ ) can be written as r = z - w, where z ∼ N (0, I) and w ∼ N (0, Σ′′ ). Then, Φ(0 | -µ′ , Σ′ ) can be approximated as

where samples {w(k) } are subject to N(µ′ , Σ′′ ) and Φ is the CDF of the standard normal distribution. Note that according to (4), we have Σ′′ = Σ′ - I = UΣU - I = U(Σ - I)U.

4 Deep Hurdle Network

In this paper we provide a deep generalization of the hurdle model, within an autoencoder framework, which we call the deep hurdle network (DHN). The DHN integrates the MVP to model the joint distribution of multiple response variables being zero or positive, and the multivariate log-normal distribution (MLND) to model the positive response variables. The MVP and MLND share the same latent features, and differences between their covariance matrices are penalized. Specifically, given a dataset D = {(x, y)(i) | i = 1, . . . , N}, where x ∈ RM is the feature data and y ∈ RL0 (or y ∈ NL0 ) is the label data, the hurdle network contains three parts (see Figure 1 for an illustration):

1. Encoder: An encoder maps the raw features x ∈ RM to latent features x′ ∈ RM′ .
2. MVP: Label data y ∈ RL0 /NL0 are translated into binary label data y′ ∈ {0, 1}L, where y′ = (y′1 , . . . , y′L) and y′j equals to 1 if yj > 0 and 0 otherwise. We then use MVP to map y′ to latent variable r ∈ RL via equations (1)-(2), where r is assumed to follow N(µ,Σ). A multilayer perceptron (MLP) is then used to model µ given x′ as input. Σ is a global parameter which is learned from random initialization and shared by all data points.
3. MLND: Let y+ ∈ RP+ (P ≤ L) be the positive part of y ∈ RL , where P is the number of positive elements of y. log(y+) = (log(y+ 1 ), . . . , log(y+ P )) is directly modeled as a multivariate normal distribution, and y+ are assumed to follow a multivariate log-normal distribution. Therefore, we have

where µ′+ and Σ′+ are the sub-parts of µ′ ∈ RL and Σ′ that correspond to y+ respectively. Note that here Σ′ is encouraged to be similar to the convariance matrix Σ of the MVP and µ′ ∈ RL is modeled with an MLP.

On the other hand, let y+ ∈ NP+(P ≤ L) be the positive part of y ∈ NL. Each y+ j is assumed to follow an univariate Poisson distribution Pr(y+ j | λ+j ) for 1 ≤ j ≤ P, where λ+j is the mean, and log(λ+) = (log(λ+1 ), . . . , log(λ+P )) is assumed to follow a multivariate normal distribution. Therefore, we have

There are several advantages of the deep hurdle network over the conventional hurdle model:

1. The encoder is forced to learn the salient features and ignore the noise and irrelevant parts of the raw features. This relieves us from selecting which salient features to use in the conventional hurdle model.
2. DHN adopts MVP and MLND to handle correlations between multiple response variables explicitly via covariance matrices, which is not considered in the conventional hurdle model.
3. The two components of the conventional hurdle model are independent, while in DHN the MVP and MLND are linked by sharing the same latent features, and penalizing the different between their covariance matrices.

4.1 End-to-End Learning for DHN

Parameters of a deep model are usually estimated by minimizing the negative log-likelihood (NLL). We develop two different learning objectives for nonnegative continuous and count data, respectively. After selecting our objective function, we can estimate the parameters of the deep hurdle network by minimizing the objective function via stochastic gradient descent (SGD).

Learning Objective for Nonnegative Continuous Data

If response variables y are nonnegative reals, then we combine equations (4), (5) and (6) to obtain the negative log-likelihood (NLL) function:

The first part of the right-hand side of equation (8) can be approximated by a set of samples {w(k) } from N(Uµ, U(ΣI)U): -log( 1K ∑Kk=1 exp(∑Lj=1 log(Φ(w(k) j )))). Such a set of samples can be obtained by first using the Cholesky decomposition, which decomposes Σ - I into CC T , and then generating independent samples v(k) ∼ N(0, I) to yield w(k) = U(µ + Cv(k) ). According to the affine transformation of the normal distribution, the samples {w(k) } are subject to N(Uµ,U(Σ - I)U). However, doing the matrix multiplications Uµ and UCv(k) are unnecessary since U = diag(2y′ - 1) ∈ {-1, 0, 1}L×L is the diagonal matrix with vector 2y′ - 1 as its diagonal. Let w′(k) = µ + Cv(k) such that {w′(k) } are subject to N(µ,Σ-I), then it is easy to show that log(Φ(w(k) j )) = y′j log(Φ(w′(k) j ))+(1-y′j ) log(1Φ(w′(k) j )). Thus, we can approximate equation (8) as follows:

where samples {w(k) } are subject to N(µ, Σ - I).

Learning Objective for Count Data

If response variables y are nonnegative integers, let y+ ∈ NP+ be the positive part of y ∈ NL0 (P ≤ L), where P is the number of positive elements of y, then each y+ j is assumed to follow a univariate Poisson distribution Pr(y+ j ; λ+j ), and log(λ+) = (log(λ+1 ), . . . , log(λ+P )) is assumed to follow a multivariate normal distribution N(µ′+, Σ′+). The likelihood of response variables y+ is given as

We then combine equations (4), (5) and (10) to obtain the NLL function as follows:

where samples {w(k) } and {log(λ+)(k) } are subject to N(µ,Σ- I) and N(µ′+,Σ′+), respectively.

Lastly, we add an L1 loss term to equations (9) and (11), respectively, to penalize the difference between Σ and Σ′ .

5 Experiments

In this section, we compare our model with the existing stateof-the-art baselines on two challenging real-world species distribution datasets concerning bird and fish populations.

5.1 Datasets

Sea bottom trawl surveys (SBTSs) are scientific surveys that collect data on the distribution and abundance of marine fishes. In each haul (sample), all catches of each species were weighted and recorded, as well as the time, location, sea surface temperature and depth. For each haul, we have additional environmental features from the Simple Ocean Data Assimilation (SODA) Dataset [Carton et al., 2018] and Global Bathymetry and Elevation Dataset , such as sea bottom temperature, rugosity, etc. Thus, we have a 17-dimensional feature vector for each haul. We study fish species abundance in the continental shelf around North America and consider SBTSs that were conducted from 1963 to 2015 [Morley et al., 2018]. We consider the top 100 most frequently observed species as the target species and there are 135, 458 hauls that caught at least one of these species. Among the 135, 458 ∗ 100 data entries, only 8.5% of them are nonzero. The distribution of the positive data is visualized in Figure 2.

eBird is a crowd-sourced bird observation dataset [Munson et al., 2012]. A record in this dataset corresponds to a checklist that an experienced bird observer uses to mark the number of birds of each species detected, as well as the time and location. Additionally, we obtain a 16-dimensional feature vector for each observation location from the National Land Cover Dataset (NLCD) [Homer et al., 2015] which describes the landscape composition with respect to 16 different land types such as water, forest, etc. We use all the bird observation checklists in North America in the last two weeks of May from 2002 to 2014. We consider the top 100 most frequently observed species as the target species and there are 39, 668 observations that contain at least one of these 100 species. Only 19.8% of these 39, 668 ∗ 100 data entries are nonzero. The distribution of the positive data is visualized in Figure 2.

5.2 EvaluationMetrics

For the task of multi-target regression, we employ and adapt two well known measures: the average correlation coefficient (ACC) 1 and root mean squared error (RMSE) [Borchani et al., 2015; Xu et al., 2019]. The ACC is defined as:

where L is the number of response variables, Ntest is the number of testing samples, y(i) and ˆ y(i) are the vectors of the actual and predicted outputs for x(i) , respectively, and ¯ y and ¯ˆ y are the vectors of averages of the actual and predicted outputs, respectively. Since the data considered here are zero-inflated, using standard RMSE might not be appropriate. Models can produce degenerate results by simply predicting a near-zero vector for each test data point. Therefore, we adapt the RMSE to the zero-inflated setting by considering zero and positive parts of the output separately. Let y(i) be the actual L-dimensional output vector of the i-th testing point, I0 (respectively, I+) be the set of indices of zero (respectively, positive) elements in y(i) , and ˆ y be the predicted L-dimensional output vector, then we define the zero-inflated RMSE (zRMSE) as:

where 0 ≤ α ≤ 1 is the relative importance of the zero part.

From the above new definition, we can see that the zero and positive parts of an output are both considered, therefore, a model cannot cheat by predicting near-zero vectors by ignoring the positive parts of the outputs. Finally, we also compare the time to learn different models.

1 The values range between -1 and 1, where -1 shows a perfect negative correlation, while a 1 shows a perfect positive correlation.

Table 1: Performance comparison between selected models on the SBTS (top) and eBird (bottom) datasets. The best scores are in bold. Except for the MLZILN/MLZIP model, each score is the average after 3 runs. α of zRMSE is set to be 0.5.

5.3 Baselines

We consider both the state-of-the-art hurdle/zero-inflated models in statistics and multi-target regression models in machine learning:

1) Hurdle/zero-inflated models. For the bird counting data, we select the state-of-the-art multi-level zero-inflated Poisson model (MLZIP) [Almasi et al., 2016] as a baseline. For the fish biomass data, we modify the MLZIP by replacing the Poisson with log-normal, which we denote as multi-level zero-inflated log-normal (MLZILN) model, as a baseline. For both the baselines, response variables are divided into random clusters. 2) Multi-target regression models. As we have discussed in Section 2, there are two kinds of methods for multi-target regression models: for problem transformation methods, we select the state-of-the-art multi-target regressor stacking (MTRS) [Spyromitros-Xioufis et al., 2016] as a baseline; and, for algorithm adaptation methods, we select the state-of-the-art multiobjective random forest (MORF) [Kocev et al., 2007], random linear target combination (RLTC) [Tsoumakas et al., 2014], multi-output support vector regression (MOSVR) [Zhu and Gao, 2018], and multi-layer multi-target regression (MMR) [Zhen et al., 2017] as baselines.

5.4 Implementations and Results

Our experiments were carried out on a computer with a 4.2 GHz quad-core Intel i7 CPU, 16 GB RAM and an NVIDIA Quadro P4000 GPU card. We use grid search to find the best hyperparameters for all models, e.g. learning rates, learning rate decay ratio, and tree depth. The encoders of DHN are parametrized by 3-layer fully connected neural networks with latent dimensionalities 512 and 256, and the MLPs of MVP and MLND are 2-layer fully connected neural networks with latent dimensionalities 256. The activation functions in neural networks are set to be ReLU. Neural network models are trained with 100 epochs and batches size 128 for the eBird dataset and 256 for the SBTS dataset. We randomly split the two datasets into three parts for training (70%), validating (15%) and testing (15%), respectively.

The results are shown in Table 1. As expected, MLZILN/MLZIP and DHN outperform other multi-target regression models in terms of zRMSE because the former models capture zero-inflation of the datasets. DHN has 12.8% and 30.9% lower errors than MLZILN and MLZIP, respectively. On the other hand, DHN has the best ACC on both datasets, which shows the benefit of employing and sharing the same covariance matrix to capture multi-entities correlations. In terms of model training time, MLZILN and MLZIP are far worse than other machine learning models, while DHN and MMR are the best. We also compare zRMSE of the methods under different α ratios in Figure 3. We can see that zero-inflated models tends to perform better for positive parts, while other nonzero-inflated models tend to underestimate the positive parts.

5.5 Ablation Studies

To learn the contributions of different components in DHN, we perform simple ablation studies. We modify the DHN architecture and rerun the experiments in the following cases:

1. We remove the encoder so that both the MVP and MLND directly use the raw features as input.
2. We remove both the encoder and the MVP so that only the MLND is used.
3. We do not penalize the difference between the covariance matrices of MVP and MLND.

The results are presented in Table 2. From the results of ablation case studies, we can observe that: (1) employing an encoder to learn latent features helps; (2) the model's performance drops significantly if we do not capture zero inflation of data; and (3) penalizing the difference between the covariance matrices of MVP and MLND helps to capture the ACC among multiple entities and boost the model's performance.

6 Conclusion

To understand the distribution of species across landscapes over time is a key problem in computational sustainability, which gives rise to challenging large-scale prediction problems since hundreds of species have to be simultaneously modeled and the survey data are usually inflated with zeros due to the absence of species for a large number of sites. We refer to this problem of jointly estimating counts or abundance for multiple entities as zero-inflated multi-target regression.

In this paper, we have proposed a novel deep model for zero-inflated multi-target regression, which is called the deep hurdle networks (DHNs). The DHN simultaneously models zero-inflated data and the correlation among the multiple response variables: we first model the joint distribution of multiple response variables as a multivariate probit model and then couple the positive outcomes with a multivariate log-normal distribution. A link between both distributions is established by penalizing the difference between their covariance matrices. We then cast the whole model as an end-to-end learning framework and provide an efficient learning algorithm for our model that can be fully implemented on GPUs. We show that our model outperforms the existing state-of-the-art baselines on two challenging real-world species distribution datasets concerning bird and fish populations. We also performed ablation studies of our models to learn the contributions of different components in the model.

Another related challenging problem in computational sustainability is to forecast how species distribution might be impacted due to the long-term effects of global climate change. In order to tackle the forecast challenge, our future works will consider using recurrent neural networks to improve our model such that it would be able to handle time series data better. On the other hand, we could also borrow ideas of advanced techniques for the multi-label prediction problem to further improve the performance of MVP that is used in our model, e.g. using latent embedding learning to match features and labels in the latent space.

Acknowledgments

This work was partially supported by National Science Foundation OIA-1936950 and CCF-1522054, and the work of Jae Hee Lee was supported by European Research Council Starting Grant 637277. We thank the Cornell Lab of Ornithology and Gulf of Maine Research Institute for providing data, resources and advice. Lastly, we thank Richard Bernstein for proofreading the paper.