# Double machine learning data integration for estimating population trends from opportunistic and structured survey data

Andrew Whetten, Diaz Renata M., Hochachka Wesley M., Stillman Andrew, Courtney Lynn Davis, Auer Tom, Matt Strimas-Mackey, Gutierrez Viviana, Howell Paige E., Kershner Eric, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Andrew Whetten, Diaz Renata M., Hochachka Wesley M., Stillman Andrew, Courtney Lynn Davis, et al.. Double machine learning data integration for estimating population trends from opportunistic and structured survey data. 2026. hal-05530740

HAL Id: hal-05530740

https://hal.science/hal-05530740v1

Preprint submitted on 27 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Double machine learning data integration for estimating population trends from opportunistic and structured survey data

Andrew B. Whetten*¹, Renata M. Diaz¹, Wesley M. Hochachka¹, Andrew Stillman¹, Courtney Lynn Davis²,

Tom Auer¹, Matt Strimas Mackey¹, Viviana Gutierrez¹, Paige E. Howell³, Eric Kershner³, Daniel Fink*¹

*Corresponding Authors: abw225@cornell.edu, daniel.fink@cornell.edu

¹ Center for Avian Population Studies, Cornell Lab of Ornithology, Cornell University, 159 Sapsucker Woods Rd, Ithaca, NY, 14850

² Department of Natural Resources and the Environment, College of Agriculture and Life Sciences, Cornell University, Fernow Hall, 111, 226 Mann Dr, Ithaca, NY 14853

³ U.S. Fish and Wildlife Service, Migratory Bird Program, 5275 Leesburg Pike, Falls Church, VA 22041.

**Abstract:**

Accurately quantifying changes in wildlife populations is essential for assessing ecosystem health and stability. Models for integrated estimation of population change (trends) can expand the inferential scope of spatially-explicit trends estimation while increasing precision and accuracy. Integration needs to account for different observation processes among surveys, sampling coverages (geographic, habitat-related, and temporal), and volumes of data. These differences lead to non-stationary contributions of surveys to spatially explicit trend estimates.

To obtain integrated trend estimates with identifiable contributions from each survey, we adapt the multi-study double machine learning (DML) model of Shyr et al. 2025. In our integrated DML trend model we address the potential issue of varying availability of data from each survey using a mixture distribution, sub-model that estimates survey coverage probability based on habitat. The mixture distribution directly weights influence of each survey on integrated estimates. We show through simulation studies that a spatially explicit mixture distribution can improve estimation and statistical power of our integrated DML trend model over independent survey-specific trend models. Using data from the North American Breeding Bird Survey and the eBird project, we show that differences in sampling coverage and data volume reflect regionally-varying strengths of the surveys, and our integrated DML trend model capitalizes on these strengths allowing each survey to be used to its maximum utility.

**Keywords:** Double Machine learning; R-learner; Causal Forests; Participatory Science Data; eBird; Breeding Bird Survey; Trend Estimation; Ecological Statistics
# 1. Introduction

Accurate metrics of population change are crucial for monitoring wildlife, guiding management strategies, and fostering public awareness about species before they go extinct (Rosenberg et al. 2019, Ricketts et al. 2005). Traditionally, structured surveys—characterized by standardized protocols, trained observers, and repeated sampling at fixed locations and times—have formed the backbone of population monitoring (Pocock et al. 2017, Marsh and Trenham 2008, Link and Sauer 1997). This approach offers significant advantages by using protocols that control for potential confounding factors that can bias estimates of population change. However, the time- and resource-intensive nature of structured surveys restricts survey coverage of various landscape compositions across a region (hereafter referred to as "sampling coverage").

Restricted sampling coverage limits the scope of possible ecological inference (Clifford et al. 2021, Guo et al. 2013, May 1999, Nichols & Williams 2006).

Opportunistic survey data, contributed voluntarily by large numbers of observers, offer a complementary data source with far greater volume and coverage. Yet, analyzing these data pose significant challenges due to complex sources of confounding variation including changes in observer effort, skill, or sensory acuity (Farmer et al. 2014); spatial biases in site selection (August et al. 2020, Shirey et al. 2021, Zhang et al. 2021); demographic variation among participants (Johnston et al. 2022); and external events such as the COVID-19 pandemic, short-term incentives, or improvements to equipment or species identification tools over time (Hochachka et al. 2021, Xue 2016). Integrating opportunistic survey data with data from structured surveys can help overcome the limitations of each data source, but successfully doing so requires: (i) mitigating bias in opportunistic datasets, (ii) reconciling differences in sampling coverage, and (iii) managing disparities in data volume.
Current approaches for integrating ecological survey data often use Bayesian hierarchical models that link survey-specific observation-process sub-models with a common ecological process model to handle confounding and heterogeneity (Zhao et al. 2024, Van Ee et al. 2024, Doser et al. 2022, Schindler et al. 2022). These models flexibly handle survey-specific biases and differences in sampling coverage, including nonrandom habitat representation – for example, when one survey targets urban areas while another excludes them – which, if ignored, can produce biased integrated estimates (Schindler et al. 2022). To address this, hierarchical models explicitly model the probability of site selection and other cross-survey differences using a small set of features (Conn et al. 2017, Clare et al. 2025, Silva et al. 2025). However, specifying and fitting multiple parametric sub-models can be practically challenging, especially with large sets of potentially informative features, nonlinear effects, and interactions. Moreover, subjective decisions to balance datasets of unequal size risk overshadowing valuable information from smaller but higher-quality structured surveys (Zipkin et al. 2021). While strategies such as subsampling or down-weighting larger datasets can be implemented, they are often ad-hoc and rely on subjective model-fitting choices. Ideally, integration methods should maximize the utility of each survey dataset while transparently tracing their influence on final population trend estimates.

Machine learning (ML) methods offer a powerful alternative to Bayesian approaches, as they can adaptively detect and model complex confounding and sampling heterogeneity without relying on explicit parametric assumptions. Double Machine Learning (DML) builds on this strength, providing a rigorous statistical framework for integrating multiple data sources while simultaneously addressing the three challenges listed above. DML enables the use of machine learning models to estimate unobservable parameters by controlling for confounding and
avoiding regularization bias (Wager and Athey 2018, Chernozhukov et al. 2018). This approach
has shown promise in estimating spatially explicit population trends from opportunistic survey
data, harnessing ML models to identify complex sources of confounding from large sets of
features while producing bias-corrected estimates of high-resolution trend estimates (Fink et al.
2023). Recent DML extensions facilitate parameter estimation across datasets (Shyr et al. 2025),
including structured (e.g. randomized controlled trials) and observational studies where there is
substantial between-study heterogeneity in data collection procedures, sources of bias, sampling
coverage, and data volume.

We developed an integrated DML trend model based on the mixture model framework of Shyr et al. (2025) to estimate spatially explicit trends in North American bird populations. This framework estimates trend parameters through a cross-survey mixture, building on the modular approach of DML to generate survey-specific estimates based on opportunistic or structured data, while introducing a mixture distribution to explicitly model between-survey heterogeneity based on the relative probability of survey coverage in a new location conditioned on features at that location. We use data from two surveys with distinct sampling protocols. First is the North American Breeding Bird Survey (BBS; Sauer et al. 2013), a long-standing, structured monitoring program with 60 years of spatially stratified, repeated-measures design. BBS data can be used to estimate coarse-scale avian trends (e.g., states or ecoregions; Sauer et al. 2013, Link and Sauer 1997) but lacks sufficient habitat sampling coverage for estimating spatially detailed (e.g. landscape-scale) trends across diverse landscapes, as all data are collected roadside. Second is the eBird project, one of the largest opportunistic bird surveys globally, comprising over 50% of observations in the Global Biodiversity Information Facility. eBird employs semi-structured survey protocols (Kelling et al. 2019) that allow observers to select the timing, location, and
method of data collection. Even though eBird has more diverse sampling coverage and data
volume, the opportunistic nature of eBird introduces complex and variable biases that challenge
the estimation of trends (Fink et al. 2023).

We demonstrate the utility of our integrated DML trend model for estimation of trends
from multiple surveys with significant differences in protocols and sampling coverage using both
simulations and a case study. To the best of our knowledge, the translation of DML methodology
for applications in the environmental sciences is limited (Jiang et al. 2025, Wei et al. 2024, Fink
et al. 2023), and this work is potentially the first effort at translating multi-study DML models to
complex environmental data with disparate protocols and sampling coverage. First, we
conducted a proof-of-concept simulation to highlight the importance of mixture distribution
specification and estimation. Second, we compared integrated trend estimates to estimates
generated from independent DML models to explore changes in mean-squared error (MSE)
performance and the bias-variance tradeoffs therein. Third, we evaluated performance across a
wide range of representative species and trend scenarios to assess expected performance via
simulation. Finally, we applied our model to BBS and eBird data for the Western Meadowlark
(*Sturnella neglecta*), to produce spatially explicit, high-resolution trend estimates that showcase
how differences in sampling coverage are leveraged within the integrated DML trend model.

**2. Methods Section:**

**2.1: Overview of the approach**

The DML trend model of Fink et al. 2023 can estimate spatially heterogeneous population trends
while accounting for confounding sources of intra- and inter-annual variation in opportunistic
survey data. Here, we extend this approach to accommodate data from multiple surveys
following Shyr et al. 2025. In our integrated DML trend model, we fit independent, survey-specific DML trend models using all available data from each study and then combine the resulting trend estimates using an additional sub-model, which we refer to as the “mixture distribution”. The mixture distribution models the probability that each survey samples any given locations based on environmental features and uses these probabilities to compute a weighted average of the survey-specific trend estimates across the study area. This approach provides a principled way to combine data from multiple studies based on their respective sampling coverages and avoids arbitrary pooling or omission of data while making transparent each survey’s contribution to the estimate.

### 2.1.1: DML Trend Models

Fink et al. (2023) introduced a DML framework for estimating population trends from a single opportunistic survey. In practice, the DML trend model consists of three predictive tasks: (1) a species-distribution model predicting average population size (i.e., conditional mean model); (2) a propensity score model predicting year of observation from observer behavior and environmental features to identify and adjust for confounding sources of inter-annual variation; and (3) a trend model that estimates average rate of population change over the study period. We adopt this DML trend model as the basis for our multi-survey extension.

Formally, we estimate the inter-annual trend for a species from one survey using

$$Y = \tau(W)T + f(X) + \epsilon, \quad (1)$$

where the response variable Y is an observed measure of population size (e.g., occurrence, $Y \in \{0,1\}$ or counts, $Y \in \{0, 1, 2, 3, ...\}$ of individuals. T is time in units of years and $\epsilon$ is a stochastic error term with $E[\epsilon | X]$. The parameter, $\tau$, is the trend to be estimated and describes
the average rate of change or amount of change, depending on whether $T$ is continuous or categorical, across the study period. The function $f$ is a non-parametric model for year-invariant variation in population size based on predictor variables $X = (X_1, ..., X_k)$, which may describe ecological conditions (e.g. habitat, climate) and/or observation processes (e.g. search effort, survey time of day).

To accommodate heterogeneity in population change, $\tau$ is allowed to vary conditionally as a function of feature vector $W = (W_1, ..., W_m)$. The feature vector $W$ can include spatially referenced information. In causal inference, this corresponds to heterogenous or conditional average treatment effect estimation; in statistics it parallels a varying coefficient model (Wager and Athey 2018, Hastie et al. 1993). In the DML trend model, $W$ can include year-invariant spatial features (e.g., a location's protected status) to localize population change, and also year-varying features (e.g., changes in habitat quality) to link trends to ecological processes of population change. When ML is used to estimate $\tau(W)$, interpretable machine learning techniques can be used to support exploration of underlying correlates of population change (Stillman et al. 2025).

Within the DML trend model, a propensity score model is introduced to identify and adjust for confounding sources of inter-annual variation when estimating $\tau$ from a single survey (Nie and Wager 2021, Chernozhukov et al., 2018, Rosenbaum and Rubin 1983). The propensity score model is written as

$$T = e(X) + \delta, \quad (2)$$

where $e$ is a non-parametric function of $X$ and $\delta$ is a stochastic error term with $E[\delta | X] = 0$. The DML approach, particularly the R-learner formulation from Nie and Wager (2021), solves
the outcome model (1) and the propensity score model (2) using plug-in estimates within a residual-on-residual regression (Robinson 1988). From (2), e(X) can be estimated directly from observed quantities. From (1), the conditional mean outcome function,

$$m(X) = E[Y | X] = \tau e(X) + f(X), \quad (3)$$

can be estimated directly from observed quantities, X and Y. Substituting the plug-in estimators, ŝ(X) and ſ(X), into (1) yields the residual-on-residual regression,

$$Y - \hat{m}(X) = \tau(W)(T - \hat{e}(X)) + \epsilon. \quad (4)$$

The residual on the left side of (4) removes the year-invariant effects of X, isolating inter-annual changes in the response while adjusting for confounding sources of intra-annual variation. The residual on the right side of (4) removes confounding due to non-random exposure to different years by regressing out the effects of X on the year variable, T. Finally, using ŝ and ſ are plugged into (4), or equivalently, used in the R-learner objective (Nie and Wager 2021), to train a model that estimates the conditional trends.

The R-learner formulation of the DML allows almost any statistical or machine learning method (e.g. penalized regression, DL, boosting) to be used for each of these three modeling tasks. Treating trend estimation as a sequence of prediction tasks makes the DML model well-suited to data with many potential confounders and complex intra- and/or inter-annual variation in the observation process. Under standard potential-outcomes assumptions (Rubin 1974), DML and R-learner estimators are unbiased (Nie and Wager 2021, Chernozhukov et al., 2018) and asymptotically normally distributed (Nie and Wager 2021). Empirically, the DML trend model has shown reduced confounding bias and strong performance for estimating spatially explicit
population trends with moderately large sets of confounding features (Fink et al. 2023, Johnston et al. 2025).

### 2.1.2: Integration across surveys using a mixture distribution

Shyr et al. (2025) extended the R-learner framework to estimate a common trend as the weighted average of multiple survey-specific estimates. Suppose we have data from K surveys that respectively track change in the same species population. Between-survey variation can arise because surveys cover different parts of the population or are subject to different sources of intra- and/or inter-annual confounding. Our goal is to estimate the common trend $\tau(W)$ at unobserved locations using all available data while controlling for sources of survey-specific heterogeneity.

Let $\tau(W)$ be the weighted mixture of K survey-specific estimates,

$$ \tau(W) = \sum_{k=1}^{K} \tau_k(W)p(k | Z), \quad (5) $$

where $\tau_k(W)$ is the trend estimate from k-th survey, and $p(k | Z)$ is the probability that a new search at a location with features Z would be conducted under the k-th survey with $\sum_{k=1}^{K} p(k | Z) = 1$. We refer to the surface representing $p(k | Z)$ spatially across a region as the mixture distribution. The purpose of the mixture distribution is to capture differences in survey coverage — particularly those likely to affect the quality (i.e., bias and variance) of the resulting trend estimates. Notably, the feature set Z used for the mixture distribution does not need to match the features sets X or W used in the DML.

Substituting (5) into the residual-on-residual regression (4) extends the DML trend model to the multi-survey setting expressed by

$$ Y - \sum_{k=1}^{K} m_k(X)p(k | Z) = \sum_{k=1}^{K} \tau_k(W)(T - e_k(X))p(k | Z) + \epsilon, \quad (6) $$
where $m_k(X) = E[Y | X, S]$ is the survey-specific conditional mean outcome model, and
$e_k(X) = E[T | X, S]$ is the protocol-specific propensity score model. Equation (6) explicitly
shows how protocol-specific trend estimates are combined to estimate an integrated trend.
Because the mixture distribution is the only cross-survey term in the expression, it fully
determines how information is shared between surveys.

To illustrate this, consider two limiting cases. For the first case, consider two surveys
($K = 2$) with identical sampling coverage across the study extent. If $P(k = 1|Z) =$
$P(k = 2|Z) = 0.5$ for all $Z$, information is shared equally everywhere, and the estimator reduces
to a pooled data analysis. For the second case, consider two surveys that cover completely
distinct regions, $Z_1$ and $Z_2$ with no overlap in sampling coverage. If $P(k = 1|Z_1) = 1$ or
$P(k = 2|Z_1) = 1$, no information can be shared between surveys, and the estimator reduces to
two full independent analyses. Between these extremes lies the more common case where
surveys partially overlap in sampling coverage with respect to $Z$ in those regions, the mixture
distribution facilitates information sharing and integration can improve estimation.

The integrated DML trend model accommodates several important sources of
heterogeneity that often complicate data integration. Because each survey has a survey-specific
DML trend model, the framework naturally permits spatially-varying trends and survey-specific
specifications of feature sets $W_k$. Differences in survey coverage, confounding, sampling
designs, and observation processes are handled through survey-specific outcome and propensity
score models, each with its own feature set. For example, opportunistic survey data may require
more detailed effort features that vary substantially across space and time, whereas structured
data may require more simplified or standardized effort features.
As an extension to the R-learner formulation of the DML trend model, the integrated DML trend model retains the same flexibility: any suitable statistical and machine learning methods can be used to estimate the survey-specific outcome models, propensity score models, heterogenous trend functions, and the mixture distribution itself. Shyr et al. (2025) show that, under standard potential-outcomes assumptions, the integrated DML trend model inherits the theoretical guarantees (i.e., asymptotically unbiased and normally distribution estimates) of the underlying DML approach.

### 2.1.3: The mixture distribution $p(k | X)$

Bayes' theorem provides a useful perspective for understanding how to estimate the mixture distribution $p(k | X)$. Using Bayes' rule, $p(k | X)$ can be written as:

$$p(k|Z) = \frac{p(Z|k)}{p(Z)} p(k), \quad (7)$$

where $p(Z | K)$ is the likelihood that a location with habitat features $Z$ would be sampled under the k-th survey; $P(Z)$ is the marginal distribution of habitat features across the study area; and $p(k)$ is the prior distribution of the probabilities for each of the $k$ surveys. A common and convenient specification is a uniform prior, which assigns equal weight to all surveys. Under this specification, $p(k|Z) \propto p(Z|k)$ and the mixture distribution is determined by differences in the likelihoods. A non-uniform prior, by contrast, encodes assumptions about the relative data volumes or importance of each survey.

In practice, mixture probabilities are usually unknown and are estimated from the data using a multinomial regression model with survey label as the response. How the training data are balanced across surveys plays an analogous role to choosing the prior $p(k)$. For example, if the training data are balanced across surveys (i.e., with an equal number of data points), the
multinomial model effectively behaves as though it were using a uniform prior by removing the
influence of unequal sample sizes and forcing the model to distinguish between surveys solely
based on differences in the likelihoods of landscape composition characterized by Z.

Conversely, if the model is trained on unbalanced data (i.e., with unequal sample sizes), the
resulting estimates reflect both the likelihoods and differing sample sizes among surveys,
analogous to specifying a prior proportional to survey data volume.

## 2.2: Application

### 2.2.1: BBS and eBird data

We applied the integrated DML trend model to structured survey data from the North American Breeding Bird Survey (BBS; Sauer et al. 2013) and opportunistic survey data from eBird to estimate spatially-explicit population trends. BBS is a long-running, structured roadside survey of more than 500 bird species across North America, conducted annually since 1966 with repeat-visits to assigned roadside routes once per year (Robbins et al. 1986; Sauer et al. 2013). eBird is a global, semi-structured participatory science project that collects opportunistic bird observations following flexible checklist protocols (Sullivan et al. 2009; Kelling et al. 2019). Participants choose when and where to survey, report effort metadata, and submit complete checklists of all species detected (Sullivan et al. 2009). The eBird project's survey design has strategically attracted an unprecedented and growing number of participants with a wide range of skills to sample large parts of the globe.

We used data from both surveys in North Dakota and South Dakota from 2014-2023 to generate data for all simulation studies and our data example. We excluded the year 2020 as the BBS did not run that year. We used data from 93 BBS routes from this region and all complete checklists from eBird. For the species-specific simulation studies, we used count data from 38
different species selected to span a range of environmental niches, observation processes, and
processes driving population change. This included species representing the three major breeding
biome groups in this region - grassland, wetland, and generalist species. Refer to Appendix 3 for
the list of species.

For each BBS stop and eBird checklist location, we summarized two sets of features,
including: (1) seven features describing search effort, observer behavior (Kelling et al. 2019),
and short-term temporal variation (e.g., day of year, time of day) to account for variation in
detection rates; and (2) 34 features describing the habitat composition and configuration in a
1.5km radius neighborhood around each survey location, including elevation, topography, land
& water cover, land use, hydrology, and road density. See Appendix 4 for a comprehensive list
of the features and their sources.

**2.3: Implementing the integrated estimator**

The multi-study R-learner of Shyr et al. (2025) can be implemented in different fashions
depending on application-based modeling considerations. For BBS and eBird data, we adopted a
fully modular approach in which all survey-specific DML components (i.e., mean outcome
model, propensity score model, and trend model) are estimated independently. This design
preserves an unambiguous attribution of each survey's contribution to integrated trends.
Implementation proceeded in three steps:

* Estimate study-specific DML components: Fit $e_1(X), \dots, e_K(X), m_1(X), \dots, m_K(X)$, and $\tau_1(W), \dots, \tau_K(W)$ separately for each survey.

* Estimate mixture probabilities: Fit $p(k = K|Z)$ using a pooled analysis of all surveys.
* Estimate heterogenous treatment effects: Use equation (6) to estimate the integrated trend.

The simplicity of this approach also results in a straightforward way to deal with challenge of unbalanced data in the data integration process (Zipkin et al. 2021) since the full sample from each survey can be used to estimate the survey-specific DML components without affecting estimates from the other surveys.

We implemented our integrated DML trend model using causal forests (Athey and Wager 2019). The propensity models for BBS and eBird data were fit using a regression forest with $X$, which contains all year-invariant landscape features (i.e. long-term means) and observer characteristics relating to search effort and observer skill (Kelling et al. 2015). The marginal models for BBS and eBird were fit using a regression forest with the same $X$. The treatment effect models for BBS and eBird were fit using a causal forest with both long-term means, and yearly deviations from the long-term means, for habitat features. The mixture distribution was fit using a multinomial random forest model using only long-term mean habitat features for $Z$. We also balanced the data from each survey when fitting the mixture distribution to isolate the effect of data volume differences.

A species' expected abundance can be defined as the product of the species' occurrence rate and the expected count of the species given occurrence, within a given area and time window (Zuur et al., 2009). To provide full insight into model performance, we evaluated the integrated DML trend model estimated separately for both trend in occurrence and trend in count given occurrence. To estimate year-to-year trends variable $T_i$ was defined as a vector of nine indicator functions, indicating whether a checklist was conducted in each of the years 2014,..., 2022, where 2023 was considered the baseline year. The corresponding trend parameter $\tau$ is the
the 9-dimensional vector describing the expected change in population (occurrence or count given occurrence) compared to 2023. We focus our simulation studies and data examples on estimating the trend τ<sub>9</sub>, the expected change between 2023 and 2014.

**2.4: Simulation Studies**

We conducted three simulation studies to evaluate the behavior and performance of the integrated trend estimator. To maintain realism while minimizing assumptions, all simulations followed the empirical data-generating approach of Fink et al. (2023), replicating many of the potentially complex observational and ecological processes captured in BBS and eBird data – including the intra- and inter-annual confounding present in the original data. Additional details on the simulation studies can be found in Appendix S1.

**2.4.1 Simulation Study 1: Proof of concept**

The goal of the first simulation study was to highlight the central role of the mixture distribution
in the construction and interpretation of the integrated trends. To do this, we simulated a
simplified data integration problem involving data from two surveys — one with structured and
one with opportunistic protocols — with contrasting coverage of a study region, and compared
the integrated estimates based on three different mixture distributions.

In this simulation, the true underlying trends in species counts, and the spatial distribution of
survey coverage, were constructed to vary across the region as a function of the proportion of
dense herbaceous croplands. In this scenario, the structured survey followed a BBS-like repeat
visit protocol in which 1000 locations were surveyed once per year over 9 years (2014-2023,
excluding 2020), by the same skilled participant conducting 3-min counts at each stop. The
opportunistic survey followed eBird-like semi-structured protocols, with 27000 observations
generated over the same time period, with no repeat visit structure by participants varying skill
level and survey effort. The three specifications of the mixture distribution we considered were:

1) “Uninformed” a uniform distribution with $p(k|Z) = 0.5$ for all grid cells, giving equal weight to each survey estimate everywhere. This simple baseline is ideal for integration when the surveys have identical sampling coverage.

2) “Low-resolution” a smooth geographic model of survey coverage estimated using a simple generalized additive model with low-rank interaction between latitude and longitude. This specification weights estimates based on regional differences in survey coverage.

3) “High-resolution” a random forest model estimating $p(k|Z)$ using fine-scale habitat features, including the dense herbaceous cropland landcover feature. This is the ideal specification resulting in a highly accurate estimate of the differences in survey coverage.

We calculated trend estimates from each survey independently and integrated estimates using each of the three mixing surfaces. For each set of estimates, filtered to the portion of the study region within the species' breeding range, we calculated RMSE, coverage (as the proportion of pixels for which an 80% CI of the estimate overlapped the true, simulated value), error (as the proportion of pixels for which an 80% CI of the estimate did not overlap zero and was in the opposite direction of the simulated value), and power (as the proportion of pixels for which the 80% CI of the estimate did not overlap zero and was in the same direction as the simulated value).

Further details about the implementation of Simulation 1 can be found in Appendix S1.

### 2.4.2 Simulation Study 2: Assessing local estimate performance
The second simulation study evaluates how integration affects local, i.e. 3km location-specific mean squared error performance. To do this, we generated 100 randomized realizations of the structured and opportunistic survey datasets with the same true but unobserved spatially varying trend. The structured survey followed a BBS-like repeat visit protocol in which the same set of locations were surveyed once per year over 9 years, by the same skilled participant conducting 3-min counts at each stop. The opportunistic survey followed a realistic eBird-like semi-structured protocol, where different locations were surveyed by different participants of varying skill expending varying levels of survey effort each year. To ensure that differences were not affected by differences in sample size, both datasets were constructed with the same sample size.

The survey-specific and integrated DML trend models were fit to each realization of the simulated data. To summarize the location-specific performance, we computed the mean squared error and the associated squared bias and variance at each 3km x 3km grid cell across the 100 realizations across the study area. Further details about the implementation of Simulation 2 can be found in Appendix S1.

### 2.4.3 Simulation Study 3: Assessing integrated trend estimation performance across species and trend scenarios

The third simulation study assessed the performance of the integrated trend estimates across a range of realistic scenarios to evaluate when and by how much integrating opportunistic data can improve trend estimation relative to using structured survey data alone. Since real-word population change is unknown, we tested integrated models against simulated data with known trends, for 38 species with diverse ecological niches, observation processes, and processes driving population change. For each species, we generated an ensemble of 40 scenarios with trends in either occurrence or count, varying in direction (positive/negative), and spatial pattern
(uniform/varying), with 5 replicate realizations per scenario. We calculated trend estimates and
performance metrics as for simulation 1. Further details about the implementation of Simulation
3 can be found in Appendix S1.

**2.5: Data Example — Spatially Explicit Integrated Population Trends for Western Meadowlark in North and South Dakota.**

We applied the integrated DML trend model to estimate spatially explicit population trends for Western Meadowlark from 2014-2023 in North and South Dakota using checklist data from both BBS and eBird. The roadside-only survey design of the BBS has known coverage deficiencies (Keller & Scallan 1999, Betts et al. 2007), whereas eBird's opportunistically sampling covers diverse regions and tends to fill many of the gaps in BBS data. Within the Dakotas, differences in survey coverage are anticipated for areas with high anthropogenic activity (i.e., urban and suburban areas, highways, and areas with high ecotourism) and areas of high avian biodiversity (e.g., permanent and temporary bodies of water). Such differences in coverage may produce disparities in survey-specific trend estimates between BBS and eBird.

For this analysis, we used 27739 BBS checklists with a detection rate of 0.66 for Western
Meadowlark, and 52.4% of checklists reporting counts of 2 or more individuals. We also used
39870 eBird checklists, with a detection rate of 0.38, and 26.7% of checklists reporting counts of
two or more individuals.

**3. Results:**

**3.1: Simulation studies**

The mixture distribution within our integrated DML trend model allows us to explicitly interpret
the contribution of data from each survey on integrated trend estimates. In Simulation Study 1,
the structured survey samples herbaceous cropland with a substantially higher probability than
the opportunistic survey and concentrates effort in the northeast portion of the study area (Figure
1). Both the low- and high-resolution mixing surfaces (Figure 1B) reflect these differences in
spatial coverage, with the high-resolution surface capturing much finer detail. The true
underlying trend (Figure 2A) varies systematically from the southwest to the northeast.

Independent estimates from each survey recover the trend only within areas of strong sampling
coverage and perform poorly elsewhere - i.e., the structured survey recovers the positive trend in
the northeast, but misses the negative trend in the southwest, and the reverse applies for the
opportunistic survey. Integration using any of the candidate mixing surfaces improves the root
mean squared error (RMSE) and confidence interval coverage compared to the single-survey
estimates. The high-resolution mixing surface produces a highly detailed trend map (Figure 2C)
and yields the estimates with the best performance in terms of RMSE, confidence interval
coverage, power, and error (Figure 2C; Appendix 2 Table S1). The correspondence between each
of the mixture distributions and the corresponding estimates highlights how the mixture
distribution can be used to explicitly interpret the contribution of data from each survey on
integrated trend estimates.

In Simulation Study 2, the integrated DML trend model showed substantial reductions in
MSE compared with the BBS-only model through most of the study area (Figure 3A). Summed
across the entire study region, the total MSE of the integrated trend estimates (325.8) was 19.7%
less than the total MSE of the BBS-only estimates (405.9) and 14.4% less than the eBird-only
estimates (380.3). To gain further investigate the performance of the integrated estimates, we
also plotted the differences between the integrated and BBS-only estimates for the squared bias
and variance (Figure 3B and 3C). The largest changes in performance were due to changes in the
squared bias (Figure 3B), though the direction of differences varied in direction across the study area with the largest reductions in squared bias occurring in the eastern third of the study area where the eBird sampling density was high. The variances of the integrated estimates are smaller than BBS-only estimates at almost all locations. However, the smallest reductions in variance occurred where the bias had the largest reductions (Figure 3C). In Appendix 2 Figure S1-S4, we include results from Simulation Study 2 based on data from four other species in the same study area and found similar results.

In Simulation Study 3, we evaluated the integrated model's ability to capture true trends of varying magnitude and structure across 38 species using simulated data constructed to replicate many of the complex changes in the observation process. Across all species and realizations, the integrated DML trend model consistently yielded improvements to RMSE and coverage relative to the BBS-only model (Figure 4A, 4B). For nearly all species, more than 60% of realizations showed improvements in RMSE and coverage, and most species had RMSE reductions in 75-95% of realizations (Figure 4C). Decreases in performance, when they occurred, were small. Appendix S3 provides example maps of true and estimated trends, mixture distributions, and performance comparisons among surveys and integrated estimates for the 38 species considered here.

## 3.2: Data Example – Western Meadowlark

Integrating BBS and eBird data produced trend estimates that reflected the complementary strengths of each survey's landscape coverage (Figure 5C). The BBS DML trend model estimated increasing western meadowlark counts in the southwest and strong declines in the east, with near-zero trends elsewhere (Figure 5A). In contrast, the eBird DML trend model indicated
widespread declines across most of the region except the southwest, where its estimates align
more closely with the BBS-specific DML trend model (Figure 5B).

In the Dakotas, survey coverage probabilities differed markedly between BBS or eBird
with respect to the included environmental features (Figure 5C and Appendix S2 Figure S5).
eBird coverage is highest in areas with water features (e.g., the Missouri River, Prairie Potholes,
etc.), along major roadways (e.g., the I-90 and I-94 interstate corridors and state highways), near
urban centers (e.g., Sioux Falls, SD, Fargo, ND, Grand Forks, ND, and Rapid City, SD), and
around national monuments and protected areas (e.g., Badlands National Park, Theodore
Roosevelt National Park, etc.; Appendix 2, Section 4). Many of these areas show strong
decreasing trends in the eBird-specific DML trend model.

The integrated trend estimates retained characteristics of both the BBS- and eBird-
specific DML estimates, drawing more heavily from eBird in landscapes where its coverage is
strongest (Figure 5). Consequently, the integrated trend map shows more pronounced declines
along major waterways, highways, and urban corridors – especially the Missouri River and
Prairie Pothole regions – where BBS coverage is limited.

**4. Discussion:**

We have shown that combining structured and opportunistic survey data using a multi-study
DML modeling approach can improve accuracy and precision of trend estimation. When survey
data sources provide complementary coverage of a study area, improvements can be substantial,
and in cases where there is increasing overlap in coverage, no harm is done through integration.
These improvements are directly related to the mixture distribution sub-model. If the mixture
distribution can effectively distinguish differences in survey sampling coverage (assuming they
exist), there is potential for the integrated trend model to outperform a survey-specific trend
model. Further, our model illustrates the potential of multi-study DML methodology for large
environmental data to strengthen inference from disparate data sources (Wang et al. 2021,
Loewinger et al. 2019, Parikh et al. 2023, Branter et al. 2024). Compared to existing models for
data integration in ecology, our work illustrates novel use of machine learning to better utilize
large, diverse sets of interacting features describing both ecological and observation processes.

Integrated trend estimation has been challenged by (i) interpreting how confounding
affects each data source separately and (ii) interpreting the influence of each data source on trend
results. For (i), existing data integration models often conflate intra- and inter-annual variation
within survey-specific observation sub-models, potentially complicating specification and
interpretation of these processes and how they affect estimates (Silva et al. 2025, Zhao et al.
2024, Conn et al. 2017, Doser et al. 2022). In contrast, the modular design of the integrated DML
trend model allows direct specification and interpretation of how within- and between-survey
variation is handled, while the use of ML methods to fit these models allows for the efficient use
of large sets of complex features data. First, the marginal model accounts for variation in
observation processes independent of years (i.e., intra-annual variation). Second, the propensity
score model describes variation in observation processes across years (i.e., inter-annual
variation). Third, the treatment effect model uses plug-in estimates from the first two models to
estimate trends conditioned on habitat features explaining spatial heterogeneity in trends across a
study region. By fitting these models independently for each survey, the approach objectively
accounts for potential differences in intra- and inter-annual variation in the observation processes
as well as spatial variation in trend estimates. Differences in survey sampling coverage are then
explicitly modeled via the mixture distribution, which transparently quantifies the relative
contributions of each data source to the integrated trend estimates. For (ii), when surveys have
disparate data volumes, many likelihood-based integration approaches can struggle, often requiring ad hoc modeling choices (Zipkin et al. 2021). Additionally, these approaches lack transparency about the contribution of individual surveys, instead relying on indirect scaling parameters, preferential sampling, or weighting strategies (Zhao et al. 2024, Silva et al. 2025, Conn et al. 2017).

Recently, there has been growing interest in assessing how well integrated trend models developed in one study system can be used to infer patterns in less-studied regions (Patterson et al. 2025). This is one line of study where the mixing surface can offer critical benefits in practice. In our integrated trend model, transfer of inference is controlled by sampling coverage, which can be explicitly modeled through the mixture distribution. We can precisely identify where the sampling coverage of one survey is not representative – and therefore not transferable – and restrict the influence of the survey in those regions. This approach parallels recent treatment effect literature on multi-study models that adaptively weight estimates based on site-level characteristics (Han et al. 2025, Han et al. 2023, Liu et al. 2023).

The Western Meadowlark data example illustrates how the mixture distribution regulates influence of BBS and eBird survey data across a large region. The mixture distribution detects areas where eBird coverage probability is higher, and the integrated model capitalizes on eBird's strong sampling coverage. For this example, differences in sampling coverage align with some differences in the estimated trend. Interestingly, the integrated trends estimated for this region align closely with the historical status and trend records for Western Meadowlark in North Dakota. The North Dakota Game and Fish Department has described Western North Dakota as a primary range for the species with generally stable population while the eastern portion of the state has been described as a secondary range experiencing declines from 1967-2012 (Dyke et al.
2015; pg. 263-264). Our results indicate that a similar trend is continuing in the Eastern and central portion of the state and only a small portion of North Dakota and South Dakota does not show evidence of a stronger declining trend.

Simulation Study 3 provides an initial view of the performance of integrated DML trend models across a suite of 38 species. In general, integration either improves or does not degrade performance relative to estimates made using BBS data alone. However, there is clear variation in performance gains among species and simulation scenarios in how much and how consistently integration improves performance. An important next direction for this work will be to formalize potential indicators of data situations in which integration is (i) merited or not, and (ii) and when it is likely to have the biggest impact. For example, low detection rates and weak count signals can make it challenging to accurately estimate trends, especially for survey-specific analyses. It is important to avoid integration in these situations, because of the potential to degrade performance. In cases where one survey has insufficient data to independently generate robust trend estimates, it may be advantageous to integrate data from both surveys when estimating the DML component models (propensity score, conditional mean model and treatment model). However, this kind of component-wise integrative modeling would inherently reduce the interpretability of survey influence on integrated trends through the mixture distribution.

Another possible strategy with the integrated DML trend model is to consider an expert-informed mixture distribution in which influence of a survey can be limited to smaller areas where data quality is known to be higher or where a stakeholder wishes to zero-out influence of a particular survey. Similarly, the mixture model framework can be applied post-hoc to integrate previously published trend estimates (along with their uncertainty estimates), similar to a meta-analysis approach to combining information across published studies.
**5. Conclusion:**

Given that spatial scale of studies remains a central challenge in ecology (May 1999), the ability of our integrated DML model to quilt together inference from multiple sources via data-driven and/or expert-informed mixture distributions represents a major contribution to statistical ecology. We have proposed and evaluated fine-scale integrated trend estimation over large geographic regions using double machine learning (DML). The integrated DML trend model presented here delivers a methodology for generating robust fine-scale trend estimates by combining spatially stratified repeat-visit survey data (e.g., BBS) with opportunistic survey data (e.g., eBird), thereby leveraging the massive volumes and spatial coverage of opportunistic data to enhance the spatial detail and geographic scope of inference derived from structured surveys.

The DML framework effectively addresses major challenges in harmonizing these diverse data sources. With the growing need to estimate continental-scale population trends, integrated trend estimation using multiple surveys holds great promise for producing stronger, more comprehensive inferences.

**USFWS Disclaimer:**

The findings and conclusions in this article are those of the author(s) and do not necessarily represent the views of the U.S. Fish and Wildlife Service.

**Data Availability Statement:**

Data and code to reproduce the analyses supporting this manuscript are available at this link: https://cornell.box.com/s/bdkrtbpsn6a39m9hwkb80vgylaztjr4z. This will be archived to Zenodo or another appropriate repository upon manuscript acceptance.

**Acknowledgements:**
We are grateful to the many thousands of participants, reviewers, and partner organizations around the world who support and contribute to eBird. We thank the eBird and engineering teams for their tireless efforts on the project. This work would not be possible without the generous support of the Cornell Lab of Ornithology's members, donors, and many staff teams.

**Funding Statement:**

This work was supported by the U.S. Fish and Wildlife Service, Grant/Award number F23AC03035-00 and by The Wolf Creek Foundation. This work used Bridges2 at Pittsburgh Supercomputing Center and Anvil at the Rosen Center for Advanced Computing at Purdue University through allocation DEB200010 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.

**Conflict of Interests Statement:**

The authors have no conflicts of interest to declare.

**References:**

1. Athey, S., & Wager, S. (2019). Estimating treatment effects with causal forests: An application. *Observational studies, 5*(2), 37-51.

2. August, T., Fox, R., Roy, D. B., & Pocock, M. J. (2020). Data-derived metrics describing the behaviour of field-based citizen scientists provide insights for project design and modelling bias. *Scientific reports, 10*(1), 11009.

3. Betts, M. G., Mitchell, D., Diamond, A. W., & Bêty, J. (2007). Uneven rates of landscape change as a source of bias in roadside wildlife surveys. *The Journal of Wildlife Management, 71*(7), 2266-2273.

4. Brantner, C. L., Nguyen, T. Q., Tang, T., Zhao, C., Hong, H., & Stuart, E. A. (2024). Comparison of methods that combine multiple randomized trials to estimate heterogeneous treatment effects. *Statistics in medicine, 43*(7), 1291-1314.
5. Bowler, D. E., Callaghan, C. T., Bhandari, N., Henle, K., Benjamin Barth, M., Koppitz, C., ... & Bonn, A. (2022). Temporal trends in the spatial bias of species occurrence records. *Ecography*, **2022**(8), e06219.

6. Burner, R. C., Kirschbaum, A., Hostetler, J. A., Ziolkowski Jr, D. J., Anich, N. M., Turek, D., ... & Niemuth, N. D. (2024). Informative priors can account for location uncertainty in stop-level analyses of the North American Breeding Bird Survey, allowing fine-scale ecological analyses. *Ornithological Applications*, **126**(4), duae041.

7. Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. *Econometrics Journal*, **21**(1):1-68.

8. Clare, J. D., Zuckerberg, B., Nunes, L. A., Strange, J., Hatfield, R., & Gratton, C. (2025). Volunteers Sample Where Endangered Bumble Bees Occur: Model-Based Triage of Preferential Sampling in Multi-Species or Integrated Distribution Models. *Diversity and Distributions*, **31**(5), e70034.

9. Clifford, S., Sheagley, G., & Piston, S. (2021). Increasing precision without altering treatment effects: Repeated measures designs in survey experiments. *American Political Science Review*, **115**(3), 1048-1065.

10. Conn, P. B., Thorson, J. T., & Johnson, D. S. (2017). Confronting preferential sampling when analysing population distributions: diagnosis and model-based triage. *Methods in Ecology and Evolution*, **8**(11), 1535-1546.

11. Dyke, Steve R., Sandra K. Johnson, and Patrick T. Isakson. 2015. North Dakota State Wildlife Action Plan. North Dakota Game and Fish Department, Bismarck, ND.

12. Doser, J. W., Leuenberger, W., Sillett, T. S., Hallworth, M. T., & Zipkin, E. F. (2022). Integrated community occupancy models: A framework to assess occurrence and biodiversity dynamics using multiple data sources. *Methods in Ecology and Evolution*, **13**(4), 919-932.

13. Farmer, R. G., Leonard, M. L., Mills Flemming, J. E., & Anderson, S. C. (2014). Observer aging and long-term avian survey data quality. *Ecology and Evolution*, **4**(12), 2563-2576.
14. Fraixedas, S., Lindén, A., Piha, M., Cabeza, M., Gregory, R., & Lehikoinen, A. (2020). A state-of-the-art review on birds as indicators of biodiversity: Advances, challenges, and future directions. *Ecological Indicators*, 118, 106728.

15. Fink, D., Johnston, A., Strimas-Mackey, M., Auer, T., Hochachka, W. M., Ligocki, S., ... & Rodewald, A. D. (2023). A double machine learning trend model for citizen science data. *Methods in Ecology and Evolution, 14*(9), 2435-2448.

16. Guo, Y., Logan, H.L., Glueck, D.H. et al. Selecting a sample size for studies with repeated measures. *BMC Med Res Methodol* **13**, 100 (2013). https://doi.org/10.1186/1471-2288-13-100

17. Han, L., Hou, J., Cho, K., Duan, R., & Cai, T. (2025). Federated adaptive causal estimation (face) of target treatment effects. *Journal of the American Statistical Association*, 1-14.

18. Han, L., Shen, Z., & Zubizarreta, J. (2023). Multiply robust federated estimation of targeted average treatment effects. *Advances in Neural Information Processing Systems*, 36, 70453-70482.

19. Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical learning. Springer Series in Statistics.

20. Hastie, T. and Tibshirani, R., 1993. Varying-coefficient models. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 55(4), pp.757-779.

21. Hochachka, W. M., Alonso, H., Gutiérrez-Expósito, C., Miller, E., & Johnston, A. (2021). Regional variation in the impacts of the COVID-19 pandemic on the quantity and quality of data collected by the project eBird. *Biological Conservation*, 254, 108974.

22. Isaac, N. J., Jarzyna, M. A., Keil, P., Dambly, L. I., Boersch-Supan, P. H., Browning, E., ... & O'Hara, R. B. (2020). Data integration for large-scale models of species distributions. *Trends in Ecology & Evolution*, 35(1), 56-67.

23. Jiang, Y., Li, L., & Xu, Y. (2025). Can digital economy improve urban ecological development? evidence based on double machine learning analysis. *Frontiers in Environmental Science*, 13, 1542363.

24. Johnston, A., Rodewald, A. D., Strimas-Mackey, M., Auer, T., Hochachka, W. M., Stillman, A. N., ... & Fink, D. (2025). North American bird declines are greatest where species are most abundant. *Science*, 388(6746),
25. Johnston, A., Hochachka, W. M., Strimas-Mackey, M. E., Ruiz Gutierrez, V., Robinson, O. J., Miller, E. T., ... & Fink, D. (2021). Analytical guidelines to increase the value of community science data: An example using eBird data to estimate species distributions. *Diversity and distributions*, 27(7), 1265-1277.

26. Kaizar, E. E. (2015). Incorporating both randomized and observational data into a single analysis. *Annual Review of Statistics and Its Application*, 2(1), 49-72.

27. Keller, C. M., & Scallan, J. T. (1999). Potential roadside biases due to habitat changes along breeding bird survey routes. *The Condor*, 101(1), 50-57.

28. Kelling, S., Johnston, A., Bonn, A., Fink, D., Ruiz-Gutierrez, V., Bonney, R., ... & Guralnick, R. (2019). Using semistructured surveys to improve citizen science data for monitoring biodiversity. *BioScience*, 69(3), 170-179.

29. Kelling, S., Johnston, A., Hochachka, W. M., Iliff, M., Fink, D., Gerbracht, J., ... & Yu, J. (2015). Can observation skills of citizen scientists be estimated using species accumulation curves? *PloS one*, 10(10), e0139600.

30. Link, W. A., & Sauer, J. R. (1997). Estimation of population trajectories from count data. *Biometrics*, 488-497.

31. Liu, Y., Levis, A. W., Normand, S. L., & Han, L. (2024). Multi-source conformal inference under distribution shift. *Proceedings of machine learning research*, 235, 31344.

32. Loewinger, G., Patil, P., Kishida, K. T., & Parmigiani, G. (2019). Multi-Survey Learning for Real-time Neurochemical Sensing in Humans using the “Survey Strap Ensemble”. *bioRxiv*, 856385.

33. Marsh, D. M., & Trenham, P. C. (2008). Current trends in plant and animal population monitoring. *Conservation Biology*, 22(3), 647-655.

34. May, R. (1999). Unanswered questions in ecology. *Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences*, 354(1392), 1951-1959.

35. Miller, D. A., Pacifici, K., Sanderlin, J. S., & Reich, B. J. (2019). The recent past and promising future for data integration methods to estimate species' distributions. *Methods in Ecology and Evolution, 10*(1), 22-37.

36. Nichols, J. D., & Williams, B. K. (2006). Monitoring for conservation. *Trends in ecology & evolution*, 21(12), 668-673.
37. Nie, X., & Wager, S. (2021). Quasi-oracle estimation of heterogeneous treatment effects. *Biometrika*, 108(2), 299-319.

38. Pacifici, K., Reich, B. J., Miller, D. A., Gardner, B., Stauffer, G., Singh, S., ... & Collazo, J. A. (2017). Integrating multiple data sources in species distribution modeling: a framework for data fusion. *Ecology*, 98(3), 840-850.

39. Parikh, H., Morucci, M., Orlandi, V., Roy, S., Rudin, C., & Volfovsky, A. (2025). A Double Machine Learning Approach for Combining Experimental and Observational Studies. *Observational Studies*, 11(3), 249-300.

40. Patterson, C. R., Zheng, X., Foster, S. D., & Helmstedt, K. J. (2025). Transferring species distribution models to understudied systems: Does data integration help?. *Methods in Ecology and Evolution*.

41. Pocock, M. J. O., Tweddle, J. C., Savage, J., Robinson, L. D., & Roy, H. E. (2017). The diversity and evolution of ecological and environmental citizen science. PLoS ONE, 12(4), e0172579. https://doi.org/10.1371/journal.pone.0172579

42. Ricketts, T.H., Dinerstein, E., Boucher, T., Brooks, T.M., Butchart, S.H., Hoffmann, M., Lamoreux, J.F., Morrison, J., Parr, M., Pilgrim, J.D., et al. (2005). Pinpointing and preventing imminent extinctions. *Proceed. Na. Acad. Sci.102*, 18497-18501.

43. Robbins, C. S., Bystrak, D., & Geissler, P. H. (1986). *The Breeding Bird Survey: its first fifteen years, 1965-1979* (No. FWSPUB157).

44. Robinson, P. M. (1988). Root-N-consistent semiparametric regression. *Econometrica: journal of the Econometric Society*, 931-954.

45. Rosenberg, K. V., Dokter, A. M., Blancher, P. J., Sauer, J. R., Smith, A. C., Smith, P. A., Stanton, J. C., Panjabi, A., Helft, L., & Parr, M. (2019). Decline of the North American avifauna. *Science*, 366(6461), 120-124.

46. Rosenbaum, P. R., & Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. *Biometrika*, 70(1), 41-55.

47. Rubin, D. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies. *Journal of educational Psychology*, 66(5), 688.
48. Sauer, J. R., Link, W. A., Fallon, J. E., Pardieck, K. L., & Ziolkowski Jr, D. J. (2013). The North American breeding bird survey 1966–2011: summary analysis and species accounts. *North American Fauna*, (79), 1-32.

49. Shirey, V., Belitz, M. W., Barve, V., & Guralnick, R. (2021). A complete inventory of North American butterfly occurrence data: narrowing data gaps, but increasing bias. *Ecography*, 44(4), 537-547.

50. Shyr, C., Ren, B., Patil, P., & Parmigiani, G. (2025). Multi-study R-learner for estimating heterogeneous treatment effects across studies using statistical machine learning. *Biostatistics*, 26(1), kxf040.

51. Schindler, A.R., Cunningham, S.A., Schafer, T.L.J. et al. Joint analysis of structured and semi-structured community science data improves precision of relative abundance but not trends in birds. *Sci Rep* 12, 20289 (2022). https://doi.org/10.1038/s41598-022-23603-0

52. Schuessler, M., Sverdrup, E., & Tibshirani, R. (2025). Statistical Learning for Heterogeneous Treatment Effects: Pretraining, Prognosis, and Prediction. *arXiv preprint* arXiv:2505.00310.

53. Silva, D., Menezes, R., Araújo, G., Machado, A., Rosa, R., Moreno, A., ... & Garrido, S. (2025). A Zero-Inflated Spatio-Temporal Model for Integrating Fishery-Dependent and Independent Data under Preferential Sampling. *arXiv preprint* arXiv:2509.09336.

54. Silverman, S. L. (2009). From randomized controlled trials to observational studies. *The American journal of medicine*, 122(2), 114-120.

55. Stillman, A. N., Davis, C. L., Dunham, K. D., Ruiz-Gutierrez, V., Rodewald, A. D., Johnston, A., ... & Fink, D. (2025). A Framework for Assessing the Habitat Correlates of Spatially Explicit Population Trends. *Diversity and Distributions*, 31(5), e70030.

56. Sun, C. C., Hurst, J. E., & Fuller, A. K. (2021). Citizen science data collection for integrated wildlife population analyses. *Frontiers in Ecology and Evolution*, 9, 682124.

57. Thomas, L. (1996). Monitoring long-term population change: why are there so many analysis methods?. *Ecology*, 77(1), 49-58.

58. Tulloch, A. I., Possingham, H. P., Joseph, L. N., Szabo, J., & Martin, T. G. (2013). Realising the full potential of citizen science monitoring programs. *Biological Conservation*, 165, 128-138.
59. Van Ee, J. J., Hagen, C. A., Pavlacky Jr, D. C., Haukos, D. A., Lawrence, A. J., Tanner, A. M., ... & Hooten, M. B. (2024). Melded Integrated Population Models. *Journal of Agricultural, Biological and Environmental Statistics*, 1-31.

60. Wager, S., & Athey, S. (2018). Estimation and inference of heterogeneous treatment effects using random forests. *Journal of the American Statistical Association*, 113(523), 1228-1242.

61. Wang, L., Chen, B., Ouyang, J., Mu, Y., Zhen, L., Yang, L., ... & Tang, L. (2025). Causal-inference machine learning reveals the drivers of China's 2022 ozone rebound. *Environmental Science and Ecotechnology*, 24, 100524.

62. Wei, X., Cheng, M., Duan, K., & Kong, X. (2024). Effects of big data on PM2.5: A study based on double machine learning. *Land*, 13(3), 327.

63. Xue, Y., Davies, I., Fink, D., Wood, C., & Gomes, C. P. (2016). Behavior identification in two-stage games for incentivizing citizen science exploration. In *International Conference on Principles and Practice of Constraint Programming* (pp. 701-717). Cham: Springer International Publishing.

64. Zhang, W., Sheldon, B. C., Grenyer, R., & Gaston, K. J. (2021). Habitat change and biased sampling influence estimation of diversity trends. *Current Biology*, 31(16), 3656-3662.

65. Zhao, Q., Latif, Q. S., Nuse, B. L., Pavlacky Jr, D. C., Kilner, C. L., Ryder, T. B., & Latimer, C. E. (2024). Integrating counts from rigorous surveys and participatory science to better understand spatiotemporal variation in population processes. *Methods in Ecology and Evolution*, 15(8), 1380-1393.

66. Zipkin, E. F., Zylstra, E. R., Wright, A. D., Saunders, S. P., Finley, A. O., Dietze, M. C., ... & Tingley, M. W. (2021). Addressing data integration challenges to link ecological processes across scales. *Frontiers in Ecology and the Environment*, 19(1), 30-38.

67. Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A., & Smith, G. M. (2009). *Mixed effects models and extensions in ecology with R* (Vol. 574, p. 574). New York: Springer.
<figure><img src="image_3.png" /><figcaption>Figure 1: Simulation Study 1 data and mixture distribution specification. (Panel A) Data from opportunistic and structured checklists locations were generated with percentage herbaceous cropland being the primary axis of variation in sampling coverage. (Left) The highest percentage cropland was generally in the upper diagonal of the study region. (Center) Opportunistic survey checklists occur less frequently in these grid cells. (Right) Structured survey checklists occur more frequently in these areas. (Panel B) We considered three specifications of the mixture distribution for data integration. The mixture model estimates the influence of each survey on integrated trend estimates based on a survey's probability of covering a grid cell with respect to landscape composition variables. (Left) A uniform mixture distribution where both surveys cover all landscapes with equal probability. (Center) A smooth mixture distribution based on Generalized Additive Model (GAM) that captures the regional-scale variation in survey coverage. (Right) A high-resolution mixture distribution based on a Random Forest classification model that uses landscape composition features, including the percentage of herbaceous cropland.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Results from Simulation Study 1. (Panel A — Left) We generated a true underlying trend for a hypothetical species where blue indicates an increasing trend and red indicates a decreasing trend. (Center) The trend estimated using exclusively opportunistic survey data is predominantly negative and (Right) the trend using exclusively structured survey data is predominately positive. (Panel B) Using the mixture distributions illustrated in Figure 1, we estimated integrated trends for a uniform (left), smooth (center), and high resolution (right) mixture distribution. The integrated DML model using high resolution mixture distribution is most effective at capturing underlying spatial structure of the true trend surface. (Panel C) Using a mixture distribution informed by probability of survey coverage with respect to landscape composition improved coverage and estimator performance of the integrated trend model.</figcaption></figure>
Wesmea difference in MSE

<figure><img src="image_5.png" /></figure>

Difference in bias^2

<figure><img src="image_6.png" /></figure>

Difference in variance

<figure><img src="image_7.png" /></figure>

Figure 3: Results from Simulation Study 2 — heterogeneity in local estimate performance. This plot shows maps of pixel-wise (3km) differences in (Left) mean squared error (MSE), (Middle) squared bias, and (Right) variance, between integrated and BBS-only trend estimates. For all maps, decreases are shown in red and indicate better performance for the integrated estimate compared to the BBS-only estimates. Summed across the entire study region, the total MSE of the integrated trend estimates (325.8) was 19.7% less than the total MSE of the BBS-only estimates (405.9) and 14.4% less than the eBird-only estimates (380.3). Although performance improves at most locations within the study area, there are more remote portions of Western North Dakota and Central South Dakota where BBS-only performance is better compared to the integrated estimates. These regions had moderately higher probability of BBS coverage compared to eBird coverage. The largest changes in performance were due to changes in the squared bias (Center), which varied in direction across the study area. The variances of the integrated estimates are smaller than BBS-only estimates at almost all locations.
A

<figure><img src="image_8.png" /></figure>

B

<figure><img src="image_9.png" /></figure>

C

<figure><img src="image_10.png" /></figure>

D

<figure><img src="image_11.png" /><figcaption>Figure 4: Results from Simulation Study 3 — The integrated double machine learning (DML) trend model generally improves model performance with respect to root mean squared error (RMSE) across a range of species and trend scenarios. We conducted simulations with respect to combinations of point count data from 38 species in our study areas, varying in trend magnitude (high/low), direction (positive/negative), and spatial patterns (constant/varying), each repeated over 5 realizations for each scenario to explore realization-to-realization variability. (Panel A and B) Across all realizations, RMSE and coverage of the true underlying trend for the integrated DML trend model were improved over the BBS DML trend model. (Panel C and D) With respect to species, the proportion of realizations where the integrated trend model improved RMSE and coverage was almost exclusively above 50%, where the average proportions of realizations with improvements was greater than 75%.</figcaption></figure>
A

BBS estimates

<figure><img src="image_12.png" /></figure>

B

eBird estimates

<figure><img src="image_13.png" /></figure>

C

Mixing surface

<figure><img src="image_14.png" /></figure>

D

Integrated estimates

<figure><img src="image_15.png" /></figure>

Figure 5: Results from the BBS and eBird Western Meadowlark integrated DML trend model for count. (Panel A) The BBS-only double machine learning trend model estimates a strong positive trend in the southwestern portion of the Dakotas and near zero trend elsewhere. (Panel B) The eBird-only double machine learning model captures portions of the same positive trend in the southwestern portion of the Dakotas, but the eBird model captures strong negative trends in pockets of this region and widespread through the Prairie pothole regions of the eastern Dakotas and along the Missouri River through the center of the Dakotas. (Panel C) The mixture distribution illustrates the probability each 3km grid cell would be surveyed by BBS, compared to eBird, based on the landscape composition within cells. The mixture distribution shows highest probability of BBS (Green) near in agricultural and grasslands. The highest probability of eBird coverage (Purple) is near water, major roadways, urban areas, and ecotourism.