Three-quarter Sibling Regression for Denoising Observational Data

1 Introduction

Observational data is increasingly important across a range of domains and may be affected by measurement error. Failure to account for measurement error may lead to incorrect inferences. For example, instrument noise in telescope data can prevent detections of exoplanet transits [Sch¨olkopf et al., 2015]; under-reporting of drug use may lead to biased public health decisions [Adams et al., 2019]; label noise in machine learning training data may lead to suboptimal models [Nettleton et al., 2010; Fr´enay and Verleysen, 2014]; and imperfect detection in ecological surveys may lead to incorrect conclusions about species populations and demographics without the proper modeling [MacKenzie et al., 2002; Hutchinson et al., 2017]. It is therefore important to develop statistical approaches to model and correct for measurement errors.

This paper is motivated by the analysis of ecological survey data. Surveys conducted by humans provide information about population sizes and dynamics for scientific understanding of animal populations and for setting conservation policies. However, humans may fail to detect a species for a range of reasons, including animal behaviors, weather, and observer skill. Making decisions to sustain animal populations requires correctly interpreting survey data with these sources of error.

∗Contact Author

There are several lines of existing work for handling measurement error. Prior work in the ecology literature has focused on explicitly modeling the detection process using latent-variable models [MacKenzie et al., 2002; Royle, 2004] . These generally assume some repetition in surveys, which helps distinguish variability in the detection process from intrinsic variability in animal counts. As citizen-science data increases in importance, another line of work has sought to model observer variability, either directly in latent variable models [Yu et al., 2010; Yu et al., 2014a; Hutchinson et al., 2017], or as separate metrics [Yu et al., 2014b; Kelling et al., 2015], for example, to help explain observerrelated variability in regression models. Latent-variable models are also used to model observation error in public health and machine learning with label noise [Fr´enay and Verleysen, 2014]. Under certain parametric assumptions, the parameters of latent-variable models for measurement error are identifiable-meaning it is possible to correctly attribute variability to the detection process as opposed to the underlying process-even without repeated surveys [Lele et al., 2012; S´olymos and Lele, 2016; Adams et al., 2019]. However, there is also vigorous debate in the same literature about the assumptions required for identifiability.

Causal modeling is an appealing alternative to parametric latent-variable models. [Sch¨olkopf et al., 2015] presented a method called "half-sibling regression" that uses causal independence assumptions to detect and remove systematic measurement error. Independence assumptions are significantly easier to reason about than specific parametric assumptions about the relationships between variables. The basic idea is to examine simultaneous measurements of multiple quantities that are known a priori to be independent. Any dependence in the measurements must be due to measurement error, which can be quantified and partially removed. Although appealing, half-sibling regression only applies when the hidden variables of interest are independent. We wish to apply similar reasoning to remove systematic noise in surveys of many species. If counts for all species are lower than expected on a given day, this may well be due to detectability and not actual changes in the populations. However, there are many common factors such as habitat and time of year that influence the true counts of different species, so these variables are not independent, and half-sibling regression does not apply.

We introduce a method called "three-quarter sibling regression" (3QS-regression) that extends half-sibling regression and can remove systematic errors in measurements of quantities that have observed common causes. Three-quarter sibling regression is derived from a causal model and has a simple interpretation when applied to species counts: the residuals of one species with respect to predictions made using the causal variables are used to correct counts of another species. We prove bounds on the ability of 3QS-regression to reduce measurement error under different assumptions. We apply 3QS-regression to a moth survey data set and show that it effectively reduces measurement error caused by moon brightness, which makes it harder to attract moths to artificial lights at survey sites.

2 Background: Half-Sibling Regression

The model motivating half-sibling regression is shown in Figure 1a. Here, Z1 and Z2 are two quantities of interest that are known a priori to be independent. In the application of [Sch¨olkopf et al., 2015], these are the brightness values of two distant stars. The variables Y1 and Y2 are measurements of Z1 and Z2 that are affected by a common noise mechanism N, for example, jitter in a telescope. In this model, any dependence between Y1 and Y2 must be introduced by the noise mechanism. The goal of half-sibling regression is to estimate one target value Z1. It does so using the estimator

which has the interpretation of subtracting from Y1 the portion of Y1 that can be predicted by Y2. In practice, Y2 may be vector valued (e.g., the brightness measurements of many other stars), and the graphical model need not follow the exact form given in Figure 1a as long as Y2⊥⊥Z1 but Y2 is not independent of Y1 given N, so that Y2 contains some "signature" of the noise process that produces Y1 from Z1.

2.1 Throwing out the Baby with the Bathwater

The key limitation of half-sibling regression for our purposes occurs when Z1 and Z2 are not independent a priori, for example, due to a common cause. In our application Z1 and Z2

will represent the counts of different species in a survey area, which are affected by common factors such as the habitat and time of year (see Fig. 2). In this setting, N represents factors that affect detectability and are shared across species. In our moth example, moon brightness is key factor, which induces a second source of dependence between moth counts. More generally, N may include transient factors such as weather, or observer attributes such as skill.

If Z1 is not independent of Z2, there are two sources of dependence between Y1 and Y2: (1) the a priori dependence induced by Z1 and Z2, and (2) the dependence introduced by the common noise mechanism N. If half-sibling regression is applied in this case, the correction term E[Y1 | Y2] will remove some of the true signal, or, as described by [Sch¨olkopf et al., 2015], it will throw out (some of) the baby with the bathwater. The variable Y2 contains information about the noise mechanism, but it is unclear whether or how this can be teased apart from the a priori dependence.

3 Three-Quarter Sibling Regression

We consider the model shown in Figure 1b. This extends halfsibling regression by adding the observed variable X, which is causal for both Z1 and Z2. We call X the process covariates. A key assumption of this model is that Z1⊥⊥Z2|X. This implies two things. First, there is no direct causal link between Z1 and Z2, which is appropriate for most species pairs in our application to ecological surveys. Second, there are no unobserved common causes of Z1 and Z2. It is up to the modeler to judge the validity of this assumption, which is standard in causal modeling. 1 In our moth survey application, time-ofyear will be the single process covariate, and it is reasonable to assume that there are no other common causes.

In this model, Y1 and Y2 are now "three-quarter" siblings: they share one parent and their unshared parents are siblings. 2 Mathematically, the key assumption is Z1⊥⊥Y2 | X. One of our results will also require N⊥⊥X. Again, although we are motivated by the particular generative model of Figure 1b, the method applies to any model that meets these assumptions. For example, the symmetry of the model and presence of Z2 is not required - Y2 can be any variable that contains some information about the noise mechanism and is conditionally independent of Z1 given X.

1 If either assumption fails it will lead again to the problem of (partially) throwing out the baby with the bathwater.

2 This term is most commonly used to describe animal kinship relationships

3.1 Estimator

The "three-quarter sibling regression" estimator or 3QSestimator is:

This is equal to Y1 minus a correction term. The quantity Y1 - E[Y1 | X] is the residual after predicting Y1 using X alone. This residual is partly due to intrinsic variability in p(Z1 | X), which we want to preserve, and partly due to the effect of noise when producing the measurement Y1. The variables X and Y2 should not be predictive of the intrinsic variability but may be predictive of the measurement noise. We subtract the portion of the residual that can be predicted using X and Y2 in order to correct Y1 towards Z1.

Next we prove some results for 3QS analogous to ones in [Sch¨olkopf et al., 2015]. Under the model of Figure 1b and an additional mild assumption, the 3QS-estimator more accurately approximates Z1 than Y1 does.

Theorem 1. Assume Z1⊥⊥Y2 | X and E[Y1 | X] = E[Z1 | X]. Then

The additional assumption that E[Z1 | X] = E[Y1 | X] is satisfied by an additive zero-mean noise model, as well other models. Note that this result does not require N⊥⊥X, as implied by Figure 1b. Before proving the theorem, we give an alternative expression for ˆ Z1 that is useful for analysis.

Lemma 1. An equivalent expression for ˆ Z1 is

Proof. From Eq. (1) we have

The last line holds because E[Y1 | X] is a deterministic function of (X, Y2). For any random variable U and deterministic function g it is the case that E[g(U) | U] = g(U).

Proof of Theorem 1. First, note that E[Y1 | X] = E[Z1 | X] implies E[Y1] = E[Z1]. Then

In Eq. (3), we subtracted EZ1 and added EY1, which are equal. In Eq. (4), we used the fact that conditional variance is no more than total variance. In Eq. (5), we used the fact that Z1⊥⊥Y2 | X. In Eq. (6), we used E[Z1 | X] = E[Y1 | X] and reordered terms.

3.2 Additive Noise Model

Now, further assume the following additive form for Y1:

Here Z1 is the "true" value and f(N) an additive error term due to N. More generally, we could have Y1 = φ(Z1) + f(N) where φ is an unknown transformation. Since we can never learn such a transformation from observations of Y1, we assume the form in Eq. (7), which amounts to reconstructing the hidden variable after transforming it to the same units as Y1.

Under the additive noise model, it is possible to quantify the error of the 3QS-estimator for reconstructing Z1.

Theorem 2. Assume Z1⊥⊥Y2 | X and N⊥⊥X. Under the additive model of Eq. (7), we have

Proof.

Eq. (8) uses the additive expansion of Y1 three times. Eq. (9) uses the facts that Z1⊥⊥Y2 | X and N⊥⊥X. Eq. (10) rearranges. Then, rearranging,

Therefore,

Theorem 2 says that it is possible to reconstruct Z1-up to a constant additive offset equal to the mean of the measurement error-with squared error equal to the conditional variance of the measurement error given the observed variables X and Y2. If the measurement error is completely determined by the observed variables, then Z1 + E[f(N)] is reconstructed exactly.

Corollary 1. If there is a function ψ such that f(N) = ψ(X, Y2), then ˆ Z1 = Z1 + E[f (N)].

Proof. In this case Var [f(N) ∣∣ X, Y2 ] = Var [ψ(X, Y2) ∣∣ X, Y2 ] = 0, which implies that ˆ Z1-(Z1+E[f (N)]) = 0.

Implementation

In practice, the conditional expectations in the 3QS-estimator are unknown, but can be replaced by regression models. The 3QS-estimator used in practice is:

where ˆE[B | A] is a regression model trained to predict B from A.

In our application, the Yi variables are symmetric and represent counts of different species. We first fit a regression model ˆE[Yi | X] using only process covariates for each species. Then, define Ri = Yi - ˆE[Yi | X] to be the residual for species i, and let R-i be the vector of residuals for all other species. We use the estimator

This has a very simple interpretation. We subtract the portion of the residual for the target species that can be predicted using the residuals of the other species. This is a special case of the 3QS-estimator where the regression model using (X, Y-i) as predictors is parameterized as a function of only the residual Y-i - ˆE[Y-i | X].

4 Experiments

In this section we experimentally demonstrate the ability of 3QS-regression to remove systematic measurement error in the presence of common causes using first synthetic data and then a moth survey data set from the Discover Life project.

4.1 Synthetic Experiments

We first conducts experiments using synthetic data, for which the true value of Z1 is known. This allows us to measure the ability of 3QS-regression to recover Z1 and evaluate its performance in settings where the auxiliary measurements Y2 = Y-1 contain a varying amount of information about the noise mechanism.

Methodology

Our synthetic data generation is patterned off the experiments of [Sch¨olkopf et al., 2015] for half-sibling regression, but extended to include process covariates. We simulated n species (indexed by i) whose occurrences Yi are determined as the following function of a process covariate X ∈ R and noise variable N ∈ R:Yi = w(i)

The coefficients w(i) X and w(i) N are drawn uniformly from [-1, 1] for each species, and control the relationship between the species occurrence and the process and noise variables. The function gi is a sigmoid function with randomized parameters (to control the inflection points, etc.). Finally ϵ ∼ N (0, σ2 ϵ ) is independent noise.

We conduct two experiments to simulate decreasing the uncertainty about f1(N) given Y-1. In the first case we set σ2 ϵ = 0 and increase the number of species, each with its own noise function fi. As the number of species n increases, we get more predictors of the error f1(N). In the second experiment we fix n = 2 and f1 = f2 (i.e the effect of noise on both species is exactly the same) while varying the noise σ2 ϵ . In either case, when the conditional variance of f1(N) given (X, Y-1) reduces, Theorem 2 predicts a more accurate reconstruction.

For these experiments we used the alternate form of the estimator given in Eq. 2, which means we fit regression models ˆE[Yi|X, Y-i] and ˆE[Yi|X] and computed ˆ Zi according to Eq. 2. The estimators could be shifted from Zi by any constant offset (cf. Theorem 2), so we centered them to have zero-mean to match Zi. We conduct our experiments with two different regression algorithms-support-vector regression (SVR) and gradient boosted regression trees, using implementations and default settings from the scikit-learn package. 3 Following [Sch¨olkopf et al., 2015], we create 20 different fixed instances and measure the mean squared error (MSE) of our denoised counts ˆ Zi against true value Zi.

Results. Figures 3a and 3b show the reconstruction error as a function of the number of species n and noise variance σ2 ϵ , respectively. The figures plot mean MSE and error across the runs. Fig. 3a shows that increasing the number of species in Y-i causes error to decrease with both regression methods. This is expected, because with more predictors that are correlated with f(N), we can learn a better model for the systematic error. Fig. 3b shows that as σ2 ϵ → 0 the error also tends towards zero. This is also expected, because, in this case, f(N) becomes a deterministic function of X, Y2. Figs. 3c and 3d compare the mean MSE of half-sibling (HS) regression against our method (3QS). HS regression performs poorly due to the common dependence of Yi variables on X.

4.2 Discover Life Moth Observations

Our second set of experiments use moth survey data from the Discover Life project 4 for studying spatio-temporal variation in moth communities. This dataset consists of counts of

3 https://pypi.org/project/scikit-learn/

4 https://www.discoverlife.org/moth

Figure 4: Residuals and moon brightness vs. time. Residuals are calculated relative to a predictive model fit across several years.

Figure 5: Pairwise scatter plots of residuals for common species and moon brightness. Inset numbers are correlations.

moths of different species collected at regular intervals at different study sites by both citizen scientists and moth experts. The protocol involves using artificial light to attract moths to a white surface and then photographing and identifying each specimen. The dataset spans eight years and has been used for analysing seasonal flights [Pickering, 2015], comparative taxonomy [Pickering et al., 2016] and other studies. The data set is unique in its temporal resolution-at some of the sites, counts have been conducted almost every day for many years.

A typical use of this type of data is to estimate smooth curves of species abundance versus time such as the ones shown in Fig. 2 (see also [Dennis et al., 2013]). The smooth curve is then used as a proxy of the actual abundance at a given site in a given year to create indices of population size and to analyze temporal patterns, such as the timing of different generations (many species occur in discrete generations each year, e.g., see H. fucosa in Fig. 2). Scientists are especially interested in how populations and timing vary across years to understand population dynamics and potential links with external factors such as climate.

Moon Brightness

A known systematic bias in moth counting studies is lunar brightness. On nights when the moon is bright, moths are less likely to be attracted to the white surface where they are photographed and counted. We show evidence of this in Figs. 4 and 5. Fig. 4 shows the residuals of moth counts with respect to a fitted model of count vs. day-of-year, together with moon brightness. The residuals track the moon phases and are anti-correlated with moon brightness. For example notice the significantly positive residual during early September which seems to rise exactly as lunar brightness reaches its nadir. Fig. 5 shows pairwise scatter plots and correlation values for the residuals of three moth species together with moon brightness. All residuals are negatively correlated with moon brightness, and positively correlated with each other, with moon brightness being one contributing factor. In this example residuals are generally more correlated with each other than with moon brightness. One explanation is that moon brightness is measured using astronomical calculations and does not correspond exactly to brightness on the ground, which depends on factors like cloud cover. A second explanation is the effect of other factors besides brightness on detectability. A third (and less desirable) explanation is that there are other causal factors that are unmeasured. The reality is likely some combination of the three.

Overall Experimental Strategy

Our hypothesis is that moon brightness is one source of systematic measurement error that can be effectively detected and removed using 3QS-regression. To test this, we will use moon brightness as an external variable that is not available to our methods, and test the correlation of moth counts with brightness before and after applying denoising techniques, as well as the accuracy of predictive models on "gold-standard" test nights when lunar brightness is minimal, so this particular source of variability is removed.

Our overall goal is to develop the best smoothed models of abundance over time for individual years (cf. Fig. 2), and we hypothesize that correcting systematic errors will help. Evaluation is challenging because the true abundance in a given year is unknown. Our approach will be to fit smooth models in individual years and test their predictive accuracy for other years. It is not our hope that the predictive accuracy is perfect - indeed, we want our model to preserve variability across years. However, another source of such variability is differently aligned moon phases, which affect detection but not true abundance. A model that corrects for such systematic errors should therefore generalize better to other years by reducing this source of variation.

Methodology

We choose moth counts from one site (Blue Heron Drive) for 2013 through 2018. Lunar brightness for the time of observation is computed with the software Pyephem 5 using the latitude and longitude of the site. We use the log-transformed counts of the most common 10 species.

We hold out one year at a time for testing and make predictions using each other year, for a total of 20 train-year / test-year pairs. When predicting one test year, we first use all four other years to fit the regression models ˆE[Yi|X] and ˆE[Ri|R-i] where Ri = Yi - ˆE[Yi|X] is the residual. In all experiments, X consists only of day-of-year. We then compute ˆ Zi for each training year using the 3QS-estimator, fit a smoothed regression model ˆE[ ˆ Zi|X] to the denoised counts for a single training year at a time, and use the model to predict on the test year. This is repeated for all train-year / testyear pairs and for all species. We compare 3QS-regression to multiple baselines, including using the original measurements Y1 to estimate Z1, as well as the half-sibling (HS) regression method of [Sch¨olkopf et al., 2015]. The pyGAM package [Serv´en and Brummitt, 2018] for generalized additive models (GAMs) is used in all regression models.

Results: Moon Brightness. We first compare the relationship between counts and moon brightness before and after 3QS-regression. It is expected that most species will have negative correlation with lunar brightness initially, but the correlation will be significantly less after denoising. The results are presented in Table 1. Typical correlations between raw counts and lunar brightness are on the order of -15% to -20%; after denosing, the magnitude of correlations decrease by an average of about 7%. The HS-regression estimator is shown for comparison, and also decorrelates the counts from lunar brightness. However if counts are correlated due to a common cause, HS-regression will also remove some of the intrinsic (non-noise) variability (the "throwing out the baby with the bathwater" problem). This can also be seen in Table 1: the overall variance of ˆ ZHS is significantly reduced relative to the raw counts. We will show in our next experiment this removed variance corresponds to true variability that should not be removed.

5 https://pypi.org/project/pyephem/

Results: PredictiveAccuracy. Fig. 6a shows predictive accuracy of several methods averaged over all species and trainyear / test-year pairs. The numbers are percent improvement in MSE relative to the baseline of a GAM fit to the noisy counts in the training year. MSE is computed only on data from the test-year with moon brightness zero. "Global" is an oracle model shown for comparison. It is trained on raw counts of four training years instead of one. By fitting one model to multiple training years, it is expected to smooth out both sources of year-to-year variability (intrinsic and moon phase) and therefore predict better on a held-out year. However, this is not our modeling goal - we want to preserve intrinsic year-to-year variability and eliminate variability due to moon phase. "MB" is a model that includes moon brightness as a feature to model detection variability. The results show that 3QS outperforms all competitors, including the MB model that has access to moon brightness as a feature. The global model performs better, as expected. The HS model is worse than the baseline due to removing intrinsic variability. Fig. 6b shows the impact of using a greater number of species for Y2. MSE is reported for the full test set as well as the moon-brightness zero set. By Theorem 2, we expect more species to reduce the conditional variance and therefore improve accuracy, which is borne out in Fig. 6b on both test-sets. Accuracy on the full test set is worse than on the moon-brightness zero test set, which is expected due to extra variability in the test set.

5 Conclusion

Building on recent work on causal models, we presented three-quarter-sibling (3QS) regression to remove systematic measurement error for species that share observed common causes. We theoretically analyzed 3QS and presented empirical evidence of its success. A future line of work would be to use simultaneous measurements to quantify uncertainty, possibly combined with Bayesian methods [Ellison, 2004]. Our method may also be useful applied to other domains, such as under-reporting of drug use [Adams et al., 2019]. A final question would be to address confounding in the presence of non-independent noise, i.e., when X is not independent of N.