Estimating koala density from incidental koala sightings in South‐East Queensland, Australia (1997–2013), using a self‐exciting spatio‐temporal point process model

Abstract The koala, Phascolarctos cinereus, is an iconic Australian wildlife species facing a rapid decline in South‐East Queensland (SEQLD). For conservation planning, the ability to estimate the size of koala populations is crucial. Systematic surveys are the most common approach to estimate koala populations but because of their cost they are often restricted to small geographic areas and are conducted infrequently. Public interest and participation in the collection of koala sighting data is increasing in popularity, but such data are generally not used for population estimation. We modeled monthly sightings of koalas reported by members of the public from 1997 to 2013 in SEQLD by developing a self‐exciting spatio‐temporal point process model. This allowed us to account for characteristics that are associated with koala presence (which vary over both space and time) while accounting for detection bias in the koala sighting process and addressing spatial clustering of observations. The density of koalas varied spatially due to the heterogeneous nature of koala habitat in SEQLD, with a mean density of 0.0019 koalas per km2 over the study period. The percentage of land areas with very low densities (0–0.0005 koalas per km2) remained similar throughout the study period representing, on average, 66% of the total study area. The approach described in this paper provides a useful starting point to allow greater use to be made of incidental koala sighting data. We propose that the model presented here could be used to combine systematic koala survey data (which is spatially restricted, but more precise) with koala sighting data (which is incidental and often biased by nature, but often collected over large geographical areas). Our approach could also be adopted for modeling the density of other wildlife species where data is collected in the same manner.


| INTRODUC TI ON
Since the time of first European settlement in Australia, koalas have faced threats from humans. Threats have included hunting (until 1927 in Queensland) (Melzer et al., 2000) and fragmentation of habitat due to urbanization. An excess of koala mortalities occurs from traumatic injuries from vehicular collisions and dog attacks or diseases such as chlamydia (Adams-Hosking et al., 2011;Lunney et al., 2002;Preece, 2007). Over the past three decades, the koala population in Australia experienced a decline of 24% (Adams-Hosking et al., 2016).
To inform koala conservation efforts, it is essential to know the geographic distribution of koala habitat and koala population densities (Adams-Hosking et al., 2016;Ellis et al., 2013). While identification of suitable koala habitat areas is relatively straightforward, the estimation of koala densities in habitat suitable areas is a difficult task (MacKenzie, 2006). Usually, systematic sampling methods are used to count koalas and estimate koala population densities. These methods typically follow a defined approach to collect data from areas that are thought to be representative of the entire geographic area of interest. These include transect and distance sampling methods (Crowther et al., 2020;Dique et al., 2004;Wilmott et al., 2019) which require skilled observers to identify and count koalas (Thomas et al., 2010). For this reason, they are labor-intensive and expensive (Kjeldsen et al., 2015) and therefore only able to be carried out over small areas.
Research has described the relationship between koala tree preference and the presence of scat (Ellis et al., 2013), and koala scat prevalence has been shown to correlate well with koala population densities (Ellis et al., 1998;Lunney et al., 2009;McAlpine et al., 2006;Phillips & Callaghan, 2011;Rhodes et al., 2008). To monitor the impact of koala conservation efforts, it is important to generate long-term datasets of koala population estimates over large geographical areas rather than estimate population counts in smaller areas at infrequent intervals (Ellis et al., 1998;Lunney et al., 2014Lunney et al., , 2016. Appropriate ecological and statistical modeling techniques can help inform decision-making in wildlife conservation if empirical data are limited (Schmolke et al., 2010) and account for biases associated with the data collection method used. A recent study in South-East Queensland (SEQLD) attempted to estimate the geographic distribution of koala populations across a wide geographical area by using a statistical modeling approach to account for multiple survey methods and multiple observers and errors in the data collection process using transect survey data collected between 1996 and 2015 (Rhodes et al., 2015).
With the advancement of communication technologies and the widespread availability of dedicated mobile applications, public participation in collecting wildlife data is increasing in popularity. For instance, members of the public were invited to collect koala sighting data as part of a program titled the "Great Koala Count" in the Australian states of New South Wales and South Australia in 2012 (Sequeira et al., 2014). The Great Koala Count has generated a large amount of incidental koala sighting data using specific guidelines for data collection in preidentified geographical areas in these states.
Attempts have been made to estimate wildlife populations using incidental sighting data alone and/or in combination with survey data (Dorazio, 2014;Sequeira et al., 2014).
In SEQLD, incidental koala sighting data have been collected since 1997, although no formal field protocols are provided to members of the public to ensure that recorded observations are valid (Dissanayake et al., 2019). While these data have been used to describe trends in sighting frequency and document spatial biases associated with sighting data (Dissanayake et al., 2019), they have not, to the best of our knowledge, been used to estimate koala population densities. It had been shown that koala sightings are spatially biased toward roads (where human activity is more frequent) and sightings tend to be more common during koala breeding seasons (Dissanayake et al., 2019(Dissanayake et al., , 2021. In this study, we use SEQLD incidental koala sighting data collected over a period of 17 years to develop a modeling approach to estimate koala density, accounting for spatio-temporal detection biases and biases arising from geographic clustering of observations.

| Study area
The geographical area of interest for this study comprised 15 local government areas (LGAs) in SEQLD (Figure 1). South-East Queensland has a higher population of koalas compared with other areas of Queensland (Dique et al., 2004) but also has a relatively high level of urban development. The point locations of koala sighting events were plotted on a map of SEQLD, indicating relatively few sightings in the north and to the far west of the candidate study area.
To develop an observation window that was tractable for modeling, we constructed a convex hull around all koala sighting locations and then dilated the convex hull by 3 km to accommodate the home range of koalas sighted on the border of the observation window (de Oliveira et al., 2014). Refinement of the study area in this way reduced the number of constituent LGAs from 15 to 11 and included parts of three remaining LGAs. The selected LGAs were in the eastern and central parts of SEQLD. The northern LGA of the Fraser Coast was excluded due to an absence of koala sighting records. This reduced the study area from 57,800 square kilometers to 30,500 square kilometers.

| Koala sighting data
Our definition for a koala sighting was a live koala observed and reported to the Queensland Department of Environment and Science (DES) by a member of the public, which we call an "incidental sighting" in the remainder of this paper. Incidental koala sightings reported in SEQLD between 1997 and 2013 (n = 14, 250 sightings) were retrieved from KoalaBASE, a database of clinical admission and sighting data of koalas, developed by the University of Queensland's School of Veterinary Science (www.koala base.com. au) and managed by DES (https://envir onment.des.qld.gov.au/). The dataset for analysis included the date of each sighting event as well as the longitude and latitude of the location where the koala was actually observed.
We assumed that if two koalas were reported in the same area in the same month, these could possibly be repeated sightings of the same animal. The term "area" in this context represents the estimated home range for a koala, which was assumed to be circular with an average (A) size of 0.35 km 2 (de Oliveira et al., 2014). We developed a selection process in which we permitted a sighted koala to be seen in the same area a maximum of 12 times per year in an attempt to remove duplicate records arising from observation of the same koala at different times (i.e., this resulted in the selection of one sighting location per koala per month). To do this, we assumed the minimum distance between two koalas was given by d = 2 × √ A ÷ , or 666 m.
Observations selected per month had a minimum distance of 666 m between closest neighbors, but observations across months did not necessarily have this minimum distance to acknowledge that home ranges can overlap over time (Ellis et al., 2009).
Subsequently, ensuring this minimum distance between two koala locations per month, the sightings for each month over the 17 years were compiled, resulting in 12 months × 17 years = 204 sighting months combined into a single dataset ( Figure 2).

| Spatio-temporal modeling of koala sighting density
The density of koala sightings (providing a proxy for koala population density) was estimated as a realization of a spatio-temporal point process model (STPP) model (Baddeley et al., 2016) using koala sighting presence-only data. The STPP is a useful statistical tool that allows one to model the spatial and temporal variation of sightings within a region and time period of interest. Within this modeling framework, two main methods can be distinguished:

F I G U R E 1
Map showing the location of the study area in South-East Queensland. The superimposed rectangular land area defines the boundaries of the area used for data analysis (water areas were excluded from the analysis). The areas marked in green are the main koala conservation areas (a) mechanistic models, where subject matter knowledge is used (i.e. our current knowledge about koala habitat and koala behavior to select relevant spatial covariates, to decide an average home range of koalas, and to describe the seasonality or "clustering" of sightings in the mating versus non-mating period) to inform the probability that a sighting will occur at a particular location and a particular time point, and (b) empirical models, where the objective is to use the observed data to inform estimates of koala population density. The approach described in this paper uses a combination of the two methods.
More specifically, we assumed that the koala density at a particular location and a particular time was dependent on the three groups of variables: (a) variables that relate to spatio-temporal detection bias; (b) spatio-temporally referenced variables that are known to be associated with koala density; and (c) variables that reflect well- Two main factors were considered in the spatio-temporal modeling of koala sighting density: spatio-temporal detection bias b (x, t) and the true koala density given by q (x, t) represents the effect of observed spatio-temporal variables, is the spatial interaction between koalas, and o is an intercept term (and where t denotes year and x spatial location).
Let ℋ t denote the past history of the process (i.e., all koala sightings) up to year t. By observing ℋ t , we condition on all the koala sightings recorded up to year tand define the intensity λof the process at a location xin year tas follows: Here, we refer to (x, t | H t ) as the koala sighting density or, in other words, the expected number of sightings per square kilometer per year. We then denote the product o × q (x, t) × r(x, t | H t ), which excludes the bias b (x, t)as the (true) koala density. Intuitively, if the bias b (x, t) = 1 for all locations and years, then (x, t | H t ) and the b (x, t)< 1 or b (x, t)> 1, these would imply underestimation and overestimation of the true koala density, respectively. The term q(x, t)that incorporates the covariates that influence koala density only and the spatio-temporal detection bias b(x, t) were modeled as a log-linear regression on covariates e(x, t) and d (x, t), respectively.
In this analysis, we used the covariate distance to primary roads as d(x, t) and used all the other covariates, which we describe more in detail in the next section as e(x, t).
The third factor corresponds to clustering of koalas that cannot be explained by d (x, t) and e (x, t) was modeled as follows:

F I G U R E 2
Diagram showing the method used to select incidental koala sightings to estimate koala population density. Circles represent an average home range area for a koala. Black and red dots represent locations of koala sightings for month 1 and month 2, respectively. The black and red lines show the closest distance between two koala sighting locations. (a) shows a selection of koala sighting locations for month 1, while (b) shows a selection of koala sightings for month 2. (c) shows the two monthly datasets combined. Using this approach, monthly sightings for each year over the 17-year study period were compiled and used to estimate koala sighting density In Equation 4, the parameter (t) regulates the strength and direction of the spatial interaction between koala sightings: 0 < (t) < 1 represents an inhibitory point process; (t) = 1 is the case of no interaction, and (t) > 1 represents the situation where koala sightings are spatially aggregated. In Equation 5, u(x) is the square root of the koala home range (in kilometers squared) at location x and is a scale parameter that regulates how quickly the spatial interaction between koalas decays as a function of distance | |x − x j | | . We let (t) vary between the mating and non-mating seasons as we ex- The likelihood function given by the above equation was maximized using a numerical optimization procedure implemented in the nlimnb function in R (version 3.5.1). Since the model described above was not available in any of the existing R packages, the algorithm used in this analysis was developed from first principles in R.
We present the model results as a graph of the predicted monthly koala population over the 17-year study period and as maps of estimated koala population density (koalas per km 2 ) for each year in SEQLD, 1997-2013. We also present a map of the estimated spatial detection bias. Although the aim of this study was to estimate spatio-temporal koala density from sighting data, we also provide predicted coefficients for the covariates with 95% confidence intervals.

| Covariates
The following data were used as explanatory variables in the model: distance to primary roads (meters), land lot density (number of lots per square kilometer), mean temperature of the hottest month per year (degrees Celsius), mean temperature of the coldest month per year (degrees Celsius), precipitation of the driest month per year (millimeters), precipitation of the wettest month per year (millimeters), mean elevation above sea level (meters), and foliage projective cover (proportion).
The variable distance to primary roads was considered to only influence spatio-temporal bias associated with koala sightings (Dissanayake et al., 2019) since it provides a measure of observer access to koala habitat. The other variables were assumed to influence observed koala sightings as they represent the documented impact of climate (temperature, precipitation), environmental factors (foliage protective cover, elevation), and indicators for the presence of humans (land lot density) on koala populations.
All of the climatic variables were obtained from online spatial databases as raster maps recorded at 1 square kilometer resolution.
The climate and elevation data were downloaded from WorldClim online database (http://www.world clim.org) for the period 1950-2000 (Hijmans et al., 2005). Foliage protective cover represents the percentage of ground area occupied by the vertical projection of foliage for 2010 and was obtained from the Queensland Government spatial catalogue (https://qldsp atial.infor mation.qld.gov.au). Primary roads and land lot density were obtained as shapefiles from the same spatial database. Distance to primary roads was calculated from its shapefile. Both, distance to primary roads and land lot density were converted to raster maps.
We standardized each covariate by subtracting the mean and dividing it by the standard deviation.
The standardized covariates are presented as raster maps in Figure S1.

| RE SULTS
The  Table S1. The coefficients of land  Figure 5. The estimated koala population density in SEQLD across the whole study period (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) is shown in Figure S2 with a scale ranging from 0 to 0.04 or more koalas per km 2 . The predicted density of koalas in the study region varied throughout the study period, with a mean density estimate of 0.0019 koalas per km 2 over the 1997-2013 period (with a maximum of 0.69 kolas per km 2 ), although estimated density in the same areas remained more or less the same throughout the study period. Based on the model presented in Equation 6, θ was >0 and was small suggesting that koalas aggregate over relatively large areas.
The predicted spatial distribution of koalas is consistent with the observed sighting data, with no koalas or very low koala densities in the western part of the study area and koala densities increasing toward the eastern coast of SEQLD, with prominent pockets of relatively high koala densities in areas with (known) good koala habitat.
The percentage of the study area in SEQLD with categorized koala sighting densities (koalas per km 2 ) for each year of the 1997-2013 study period are shown in Table 1. The percentage of land areas with very low sighting densities (0-0.0005 koalas per km 2 ) remained similar throughout the study period representing a mean (SD) of 66% (0.06) of the total study area (Table 1)

| D ISCUSS I ON
We present here the results of spatio-temporal point process model where koala population density was estimated from citizen science koala sighting data by adjusting for spatio-temporal detection bias.
The mean estimated density of koalas in the study region over the 1997-2013 period was 0.0019 koalas per km 2 (with up to 0.69 koalas per km 2 ). There was marked variation in koala population density largely due to the heterogeneous nature of koala habitat across the study area. Rhodes et al. (2015) reported koala densities ranging from 0.00001 to 0.11 per km 2 in coastal regions of SEQLD, with an average of 0.0004 koalas per km 2 . Our study shows similar distribution patterns to the study by Rhodes et al. (2015), with high densities along the east coast area and two very high-density spots in the central east coast area and low densities in the western part of the study area (Dissanayake et al., 2021). However, the model developed by Rhodes et al. (2015) used data collected through multiple systematic surveys, which were implemented in small areas and did not predict koala populations across large geographic areas due to uncertainties associated with extrapolation. In fact, extrapolating koala densities from statistical models for large geographical areas is questionable as koala habitat is not continuously distributed. Ideally, to avoid this F I G U R E 5 Raster map showing estimated koala population densities (koalas per km 2 ) across South-East Queensland between 1997 and 2006. Estimates were derived from a spatio-temporal point process model using koala sighting data (n = 6,580) problem, densities should be predicted to strata of different habitat types (Dique et al., 2004), but this was beyond the scope of the methodology presented in the Rhodes et al. (2015) paper.
Actual koala numbers are difficult to estimate. In 2010, the Department of Environmental Heritage and Protection (DEHP) predicted that the Queensland koala's population was between 157,000 and 177,000 individuals, while the Threatened Species Scientific Committee of Australia estimated Queensland's koala population to be approximately 167,000 individuals in 2010, representing a 43% decline from 1990 (Rhodes et al., 2015). Another study esti- We identified a strong clustering of koalas in locations in and around the Moreton Bay and Redlands areas which is similar to the high-density areas identified by Rhodes et al. (2015) using systematic field survey data. Our model identified low densities of koalas in the western part of SEQLD whereas Rhodes et al. (2015) predicted higher densities there, although this was probably due to the uncertainty associated with the model estimates for this region.
Importantly, we were able to estimate koala population density over time and space while incorporating a range of covariates expected to be associated with observed sighting densities or spatio-temporal detection bias. For example, distance to primary roads was considered to be an explanatory variable predominately influencing spatio-temporal detection bias, while foliage protective cover was influencing an observer's ability to sight a koala and therefore impacting on observed sighting density. However, the contribution of explanatory variables to the two different components of the model cannot be quantified as these components were included as additive factors on a log-scale.
Considering that covariate coefficients with a negative sign decrease reported sightings, our model indicates that larger distances to primary roads, denser foliage, higher altitude and increased precipitation, results in decreased koala population density. In contrast, increased lot density and warmer temperatures in the colder months of the year were associated with increases in estimated koala population densities.
Uncertainties in estimated koala densities can be further reduced whether additional data are collected at the time of each sighting event. These data could then be used to estimate and remove the effect of the bias on population density estimates.
We have presented a map showing the geographic distribution of spatial detection bias ( Figure S3). Due to the formulation of our TA B L E 1 Percentage of the South-East Queensland study area, 1997-2013, belonging   In our study, no observer-related variables were collected at the time koalas were sighted. It has been shown that the probability of detection of a koala by an observer varies with their previous experience: An experienced observer can have a detection rate of around 70%, while an inexperienced observer might have a detection rate in the order of 30% (Corcoran et al., 2019). As a result, many koalas may go undetected simply because of the lack of observer experience.
The situation is somewhat different in systematically conducted field surveys carried out by trained individuals (Rhodes et al., 2015). Thus, incidental sightings reported by members of the public represent a biased sample of the koala population at any given time, but the collection of data on the experience of observers at the time of the sighting could provide valuable information to correct for this bias.
The frequency of koala sightings varies between seasons of the year. Such seasonal variations might be due to more frequent dispersal of koalas during breeding periods but also due to better visibility of animals and weather conditions that are more favorable for people to go outdoors and spot koalas. Interestingly, the results of our model indicated that the clustering of koalas is not prominently different between the mating ( 1 = 2.006) and nonmating seasons of koalas ( 2 = 2.029). This might be explained by koalas being solitary animals, and although they travel over larger distances in the breeding season, their greater mobility might not necessarily be associated with clustering of animals.
We included an average home range of koalas in our model, as we did not have detailed koala home range information for different parts of our study area. We acknowledge that koala home ranges are not uniform, and even within the Redlands Local Government Area, koala home ranges of koalas vary between 0.05 and 0.55 km 2 (de Oliveira et al., 2014). The precisions of koala densities could be improved if home ranges appropriate for each habitat type were included in the model.
It has been predicted that drier and warmer climatic conditions have an undesirable impact on koala habitat and thereby negatively impact on koala population densities (Adams-Hosking et al., 2011).
Unfortunately, our study was constrained by the nonavailability of some temporally varying explanatory variables (e.g., foliage projective cover and land lot density). Therefore, for consistency, we used a single value for each explanatory variable across all years, and, as a result, the temporal effect of explanatory variables such as the impact of temperature changes over time on koala densities could not be quantified.
In conclusion, we developed a statistical model that addressed the spatio-temporal bias associated with observed koala sightings and provided long-term koala density estimates for one of the largest koala populations of Australia over a 17-year period. Such estimates of koala population size are required for koala population viability analyses, epidemiological models, translocation programs, and other koala management programs. In future, the modeling approach presented in this paper could be used to combine systematic (spatially restricted) survey data with koala sightings data, allowing koala population densities to be estimated over larger geographical areas. The systematic survey data could also provide the means by which to validate sightings data model outputs. Repeated systematic survey data in the same location could further inform the occupancy which can be used in an extended model. Our approach could also be adopted for modeling the density of other wildlife species where incidental sightings might be recorded.

ACK N OWLED G M ENTS
We thank the Queensland Department of Environment and Science for providing access to the koala sighting data as part of the KoalaBASE project. We also thank the Moggill Koala Rehabilitation Centre staff at Brisbane for providing insights into the data collection process.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The incidental koala sighting dataset analyzed in this study was obtained from KoalaBASE (www.koala base.com.au).