Model

This section describes the zonal travel cost method, which begins with a demand model of the form:

(1) $$trip{s_{ijt}}/{N_i} = \;{e^{\rho t{c_{ij}} + {\delta _{jt}} + {\alpha _j} + {\tau _t}}}$$

where $trip{s_{ijt}}$ is the number of trips to site j from residential zone i at time t, ${N_i}$ is the population of birders living in i, and $t{c_{ij}}$ is the travel cost between i and j (Willis and Garrod Reference Willis and Garrod1991). In our application, ${\delta _{jt}}$ describes the effect of time-varying site quality, ${\alpha _j}$ is a site fixed effect, and ${\tau _t}$ is a time fixed effect. These fixed effects will control for site attributes or time-varying shocks not explicitly included in the model, such as the presence of a lake, parking, or weather. Our application pools trips across several sites and characterizes time-varying site quality using measures of species richness and crane abundance.Footnote ¹

We modify eq. (1) so that quantity demanded is trips per zone rather than trips per zonal resident:

(2) $$trip{s_{ijt}} = \;{e^{ln \left( {{N_i}} \right) + \rho t{c_{ij}} + {\delta _{jt}} + {\alpha _j} + {\tau _t}}}$$

which we estimate using Poisson regression. This means that we can interpret the parameters in terms of changes in trips rather than trips per resident. Note that eq. (2) takes the form of a count data model, which is similar to traditional single-site recreation demand models (Parsons Reference Parsons, Champ, Boyle and Brown2003). We briefly discuss this similarity further below. Although it is common in the recreation demand literature to estimate count data models using negative binomial regression to avoid problems with “overdispersion” (Holmes and Englin Reference Holmes and Englin2010), this is not necessary. Trips do not need to be Poisson distributed for the regression to produce consistent parameters, which are robust to distributional misspecification. Overdispersion can, however, severely bias Poisson standard errors. We correct this bias by reporting “robust” standard errors using the sandwich estimator recommended by Wooldridge (Reference Wooldridge1999).Footnote ²

Our approach identifies the effect of site attributes in ${\delta _{jt}}$ because we estimate eq. (2) using data pooled across months and sites. The zonal method also implicitly pools trip data across individuals coming from the same residential location. In contrast, traditional single-site demand models use data on individuals’ seasonal trips. This allows researchers to identify the effect of travel cost and individual characteristics on demand but not site attributes. So, it is not possible to value quality changes with traditional single-site demand models. However, researchers can value attributes by pooling trips over time, when quality changes, or by stacking visitation and attribute data from multiple sites (Parsons Reference Parsons, Champ, Boyle and Brown2003). Our approach does both.

We estimate several specifications for ${\delta _{jt}}$ . The first is

(3) $${\delta _{jt}} = {\beta _1}richnes{s_{jt}} + {\beta _2}anycrane{s_{jt}}$$

where $richnes{s_{jt}}$ is the count of species present and $anycrane{s_{jt}}$ is an indicator that equals one if sandhill cranes are among those present and zero otherwise. We follow Kolstoe and Cameron (Reference Kolstoe and Cameron2017) and Jayalath et al. (Reference Jayalath, Lloyd-Smith and Becker2023) by including the number of species as a site attribute.Footnote ³ The parameter ${\beta _1}$ measures the average effect of an additional species on trips per zone, while the parameter ${\beta _2}$ measures the contribution of cranes relative to other species. A positive ${\beta _1}$ would indicate that trips increase with the number of species, and a positive ${\beta _2}$ would indicate an even larger increase occurs if cranes are present. Put differently, ${\beta _2}\gt0$ implies that cranes have a larger positive effect than the presence of another species, on average, while ${\beta _2} \lt0$ means that cranes have a smaller effect than another species.Footnote ⁴

The next specification replaces the $anycrane{s_{jt}}$ indicator with a measure of abundance,

(4) $${\delta _{jt}} = {\beta _1}richnes{s_{jt}} + {\beta _2}crane{s_{jt}}.$$

where $crane{s_{jt}}$ is the population of cranes at j at time t. Crane abundance could affect trips if birders prefer sites with more rather than fewer cranes, in which case ${\beta _2}\gt0$ . Of course, by modeling this preference linearly, eq. (4) implicitly assumes that abundance has a constant marginal effect. The third specification allows for nonlinear effects by separating abundance into four bins,

(5) $${\delta _{jt}} = {\beta _1}richnes{s_{jt}} + {\beta _2}{{\rm{\Phi }}_{25jt}} + {\beta _3}{{\rm{\Phi }}_{50jt}} + {\beta _4}{{\rm{\Phi }}_{75jt}} + {\beta _5}{{\rm{\Phi }}_{100jt}}.$$

where ${{\rm{\Phi }}_{pjt}}$ is an indicator for population in quartile p, conditional on the presence of cranes. Parameters ${\beta _2}$ through ${\beta _5}$ can shed light on the effect of cranes at different population levels. For example, ${\beta _5} = {\beta _4} = {\beta _3}\gt {\beta _2}\gt0$ implies that the effect of abundance saturates when the population surpasses the median population at sites.

In the application that follows, we estimate equations (2)–(5) twice. First, we use measures of richness and crane abundance generated from eBird. Second, we replace the abundance variable with data from the Indiana DNR. The goal of eBird is to measure populations that organizations do not have the resources to track professionally, but using volunteers could lead to bias. Comparing DNR and eBird crane counts can shed light on the accuracy of transferring models and welfare estimates based on one type of data to another.

We can use the parameters to calculate consumer surplus for site access and quality changes. The formula for access value, expressed in terms of willingness to pay (WTP), is

(6) $$WTP = - {{1} \over {\rho }}$$

which is denominated in dollars per trip (Parsons Reference Parsons, Champ, Boyle and Brown2003). In our application below, we also estimate the value of a one-unit increase in species, $WT{P_{ijt}} = - {{1} \over {\rho }}\left[ {{e^{{\beta _1}}} - 1} \right]trip{s_{ijt}}$ , and a 10% increase in the number of species, $WT{P_{ijt}} = - {{1} \over {\rho }}\left[ {{e^{{\beta _1}richnes{s_{ijt}} \times 0.1}} - 1} \right]trip{s_{ijt}}$ , where the terms in brackets measure the proportional change in trips in the location and time period that experienced the change. Denominated per trip, the value of a one-unit increase is

(7) $$WTP = - {{1} \over {\rho }}\left[ {{e^{{\beta _1}}} - 1} \right],$$

and for a 10% increase,

(8) $$WTP = \mathop \sum \nolimits_i \mathop \sum \nolimits_j \mathop \sum \nolimits_t {{1} \over {\rho }}\left[ {{e^{{\beta _1}richnes{s_{ijt}} \times 0.1}} - 1} \right]/{\rm{\Psi }},$$

where ${\rm{\Psi }}$ is the number of observations. Additionally, we use the form of eq. (7) and eq. (8) to value unit and percentage changes in cranes when the effect of abundance is measured as in eq. (4).

Study area

We apply the demand model to three areas in Indiana that experience large increases in migratory species in the fall and spring. One of these species, the sandhill crane, is popular with birders because it travels in large flocks and has a loud, distinct call. Like the endangered whooping crane (Grus americana), sandhill cranes are large, standing about three feet tall with a wingspan of six feet, although unlike their biological cousin, they are not at risk of extinction. Public locations in Indiana with ideal habitat for cranes include the Jasper-Pulaski, Goose Pond, and Muscatatuck wildlife areas (A. Phelps, DNR, personal communication). These are also the properties where the Indiana DNR operates crane counts, and thus where we apply the demand model. The remainder of this section briefly describes these areas.

Jasper-Pulaski is a state fish and wildlife area that protects about 8200 acres of habitat in northwest Indiana. Developed in the 1930s as a hunting and fishing preserve, large parts of the wildlife area are now managed by the Indiana DNR as waterfowl habitat. The site includes a crane observation platform and viewing scopes. DNR biologists have recorded migration events that can bring as many as 30,000 cranes to the area on some days. The media attention this generates makes this one of the best-known birdwatching locations in the state, particularly during the fall (Greene Reference Greene2020; Zorn Reference Zorn2021).

Goose Pond is a state fish and wildlife area that includes about 9000 acres of prairie and marshland in southwest Indiana. After purchasing the property in 2005, the DNR restored a wetland complex that had originally been drained for farming. In 2016, the DNR opened a visitor center with an observation deck. The area provides habitat for a variety of wildlife, with over 260 different types of birds documented. This makes it popular with birdwatchers year-round (A. Phelps, DNR, personal communication).

Muscatatuck National Wildlife Refuge covers 7700 acres in south central Indiana. Unlike the other two sites, Muscatatuck is a federally protected area. The US Fish & Wildlife Service manages the site primarily as migratory waterfowl habitat and for wildlife viewing. Visitors have access to hiking trails and a driving loop. Similar to Goose Pond, the area provides ideal habitat for a large number of species that attracts birdwatchers throughout the year, while the number of sandhill cranes during the migration season is comparable to Jasper-Pulaski (Freedman Reference Freedman2023).

Although located in different areas, Jasper-Pulaski, Goose Pond and Muscatatuck are similar in terms of size and habitat, which consists of a mix of wetland, prairie, and bottomland forest. This similarity makes them potential substitutes for birdwatchers and good candidates for a pooled site demand model. To be clear, there are many other wildlife viewing areas in Indiana popular with birders. One concern is that excluding potential substitutes from the model could lead to bias. However, in Appendix A, we show that expanding the group of sites does not qualitatively affect the estimates. This suggests that our conclusions are robust to a larger number of choice alternatives. In fact, as we discuss below, WTP for site access and species richness from the model compare favorably with estimates from existing, more complex RUM models of birdwatching in the literature.

Data

The primary data set comes from eBird, which tracks bird observations recorded through a mobile application. After an observation or “sampling event,” volunteers submit a checklist that indicates whether birding was the primary purpose and the birds they identified by sight or sound. We downloaded the 2022 eBird database in two files: eBird data (EBD) and sampling event data (SED). The EBD records the species data from sampling events, including common and scientific name, number observed, location details, travel protocol, date, an observed identifier, a group identifier, and a sampling event ID code. The location data specifies state, county, and an eBird locality. Localities are geocoded and classified as either personal or hotspot. Personal localities, which are observer-specific, can range from an apartment balcony to the side of a county road. Hotspot localities are public areas frequently used by birders. Travel protocol records whether a volunteer made an observation while traveling, stationary, or incidentally. We used this data to drop incidental trips.Footnote ⁵ The SED contains the same variables as the EBD except that it is organized by checklist rather than species, which we used to zero-fill crane counts, as described below.

We used hotspots, the three wildlife areas and county boundaries to define sites in the model. We identified five hotspots in Jasper-Pulaski, eight in Goose Pond, and sixteen in Muscatatuck. We followed eBird’s recommendation to group hotspots in the same vicinity by aggregating those in the same wildlife area and county. This is because wildlife areas can have several hotspots with similar coordinates (e.g. different trails). Our approach ensures that localities in the same vicinity are grouped to a common site. This produced two sites for Jasper-Pulaski, one for Goose Pond, and three for Muscatatuck.

Next, we measured trips by counting the number of sampling events in the SED. To separate trips by residential zone, we assigned each event to a ZIP code tabulation area in three steps nesting an iterative procedure. The first step linked a sampling event to the ID code assigned to each eBirder. In the second step, we searched through sampling events for personal localities the eBirder named “Home” or included “home,” “residence,” etc. in the name. We used the latitude-longitude for this locality as the eBirder’s residential location. For the remainder, though, we had to interpolate home locations, using either the coordinate of an eBirder’s modal personal locality or, for anyone who never submitted a personal locality, the average latitude-longitude of visited hotspots.Footnote ⁶ This is necessary to avoid skewing the sample toward birders willing to identify their home, which could exacerbate bias. Of course, the interpolation itself could lead to bias, and in Appendix B we present estimates after applying a more conservative interpolation procedure. Third, we assigned each sampling event to the corresponding observer’s home location and hence ZIP code. ZIP codes not assigned as a residential location for any Indiana eBirder were dropped from the model. Finally, we summed trips by ZIP code and site in each biweekly period.

An important concern, as noted by Kolstoe et al. (Reference Kolstoe, Cameron and Wilsey2018), is that eBird volunteers self-report and may not be representative of the population overall. Recent research appears to validate this concern. Cameron and Kolstoe (Reference Cameron and Kolstoe2022) find that a representative sample were less willing to consider traveling to more distant sites compared with eBird volunteers. Similarly, Rosenblatt et al. (Reference Rosenblatt, Dayer, Duberstein, Phillips, Harshaw, Fulton and Wood2022) conclude that eBird volunteers are more avid than the average birder. Volunteers also have a higher average income relative to the general population (Grade et al. Reference Grade, Chan, Gajbhiye, Perkins and Warren2022), although birders on the whole tend to have above-average incomes (Carver Reference Carver2019). We address this concern using weights. First, we correct for differences in representation by geography by weighting the data by the number of households relative to the number of eBirders in a ZIP code. This stratifies the sample so that it matches the population in each zone. Second, we use the income distribution from a nationally representative sample of birders (Carver Reference Carver2019) to construct sample selection weights that account for response differences by income. Finally, we use statistics on the eBird volunteering rate as a function of birdwatching distance (Cameron and Kolstoe Reference Cameron and Kolstoe2022) to construct zonal-site weights that limit the share of individuals in a zone who would travel to more distant locations.Footnote ⁷

Next, we calculated species richness and crane abundance. After pooling the hotspot data by site, we counted the number of unique species observed in the EBD in each biweekly period, which includes the days up to the 15^th or after the 15^th of each month. We then calculated the abundance of sandhill cranes by merging the EBD and SED data sets. Most volunteers who reported seeing cranes counted the number observed. We assume an observer saw no cranes if their checklist appeared in the SED but they reported no cranes in the EBD. Thus, if a sampling event occurred at one of the hotspots (in the SED) in our study area but did not count any sandhill cranes (in the EBD), we assigned it a count of zero.Footnote ⁸ In hotspots where volunteers reported different numbers at the same time, we used the median. We then measured the eBird crane population by summing the (median observed) number of cranes across hotspots at a site. We explore alternative methods of measuring abundance in the discussion section below.

We filled in missing species richness and crane abundance data using imputation. While the study area includes hundreds of sampling events, occasionally a site had no sampling events in a biweek. We imputed these data by applying Poisson regression to the model ${s_{jt}} = {\rm{exp}}\left( {{\lambda _j} + {\mu _t}} \right)$ , where ${s_{jt}}$ is either $richnes{s_{jt}}$ or $crane{s_{jt}}$ , ${\lambda _j}$ is a county effect, and ${\mu _t}$ is a time effect. To leverage observations from other hotspots in the same counties as our study sites, which could correlate with species richness and crane abundance, we estimated the model using data on all j in Indiana, rather than just the sites in the model.

We aggregated the DNR crane counts to the same biweekly periods as the eBird data. Property staff perform these counts weekly during the migration season, which runs from the last week of August to the last week of January. The counts are made when cranes leave their roost in the morning. After estimating the number remaining in the roost, the counts are aggregated across staff. We did not impute the DNR’s missing values as we did for eBird because the DNR counts appeared complete during the migration season. To address any concerns that the timing of this season could skew the comparison between eBird and DNR data, one of our robustness checks restricts all of the eBird observations to the migration season.

We calculated travel cost using latitude-longitudes, fuel and mileage-related depreciation costs, and an estimate of the value of travel time. We calculated travel distances $dis{t_{ij}}$ and travel times $tim{e_{ij}}$ between the centroid of each ZIP code i and the average latitude-longitude of site j’s hotspots. We used AAA’s 2022 Your Driving Cost report to calculate the cost per mile of driving, by subtracting the average cost of driving 10,000 miles from the cost of driving 15,000 miles, and then dividing by 5,000. We then estimated the value of travel time by dividing the median household income $incom{e_i}$ in a ZIP code by 2000 multiplied by one-third, which assumes working 2000 hours in a year and that travel time is valued at one-third the wage rate (Parsons Reference Parsons, Champ, Boyle and Brown2003). Finally, we calculated travel cost using the formula:

$$t{c_{ij}}=2\left[ {dis{t_{ij}} \times 0.2718 + {{1} \over {3}} \times tim{e_{ij}} \times \left( {\frac{{incom{e_i}}}\over{{2000}}} \right)} \right]$$

In Appendix C, we present results after setting the value of travel time equal to one-half the wage rate.

Table 1 summarizes the data. Between six sites, 24 biweekly periods, and 540 ZIP codes, there are 77,760 observations. The average number of trips is 0.017; this is less than one because there are zero trips in most ZIP code-site-biweek combinations. The average travel cost is $117 across all observations, or $75 after weighting on trips. Using eBird, the average number of species is 43 and the average crane count is 206. The average eBird crane count is substantially less than the maximum, which exceeds 5,000. The average DNR count is 6,570 cranes and the maximum is 31,536.

Table 1.

Summary statistics of variables used in the pooled site demand model

The number of observations is 77,760 except that the DNR cranes data has 35,100 observations.

Before turning to the results, we should clarify that a trip in our data is taken primarily for the purpose of birding but could include multi-destination trips. Lupi et al. (Reference Lupi, Phaneuf and von Haefen2020) recommend excluding multi-destination trips. Unfortunately, eBird checklists do not indicate whether a single trip covered multiple sites. However, we can infer multi-destination trips by identifying eBirders who submitted more than one checklist on the same day. Twenty-three percent of trips could be multi-destination following this approach. Appendix D presents estimates after dropping these trips, which we find does not qualitatively affect the results.Footnote ⁹