Optimization of spatial scale, but not functional shape, affects the performance of habitat suitability models: a case study of tigers (Panthera tigris) in Thailand

Species habitat suitability models rarely incorporate multiple spatial scales or functional shapes of a species’ response to covariates. Optimizing models for these factors may produce more robust, reliable, and informative habitat suitability models, which can be beneficial for the conservation of rare and endangered species, such as tigers (Panthera tigris). We provide the first formal assessment of the relative impacts of scale-optimization and shape-optimization on model performance and habitat suitability predictions. We explored how optimization influences conclusions regarding habitat selection and mapped probability of occurrence. We collated environmental variables expected to affect tiger occurrence, calculating focal statistics and landscape metrics at spatial scales ranging from 250 m to 16 km. We then constructed a set of presence–absence generalized linear models including: (1) single-scale optimized models (SSO); (2) a multi-scale optimized model (MSO); (3) single-scale shape-optimized models (SSSO) and (4) a multi-scale- and shape-optimized model (MSSO). We compared performance and resulting prediction maps for top performing models. The SSO (16 km), SSSO (16 km), MSO, and MSSO models performed equally well (AUC > 0.9). However, these differed substantially in prediction and mapped habitat suitability, leading to different ecological understanding and potentially divergent conservation recommendations. Habitat selection was highly scale-dependent and the strongest relationships with environmental variables were at the broadest scales analysed. Modelling approach had a substantial influence in variable importance among top models. Our results suggest that optimization of the scale of resource selection is crucial in modelling tiger habitat selection. However, in this analysis, shape-optimization did not improve model performance.

Objectives We provide the first formal assessment of the relative impacts of scale-optimization and shapeoptimization on model performance and habitat suitability predictions. We explored how optimization influences conclusions regarding habitat selection and mapped probability of occurrence. Methods We collated environmental variables expected to affect tiger occurrence, calculating focal statistics and landscape metrics at spatial scales ranging from 250 m to 16 km. We then constructed a set of presence-absence generalized linear models including: (1) single-scale optimized models (SSO); (2) a multi-scale optimized model (MSO); (3) singlescale shape-optimized models (SSSO) and (4) a multiscale-and shape-optimized model (MSSO). We compared performance and resulting prediction maps for top performing models. Results The SSO (16 km), SSSO (16 km), MSO, and MSSO models performed equally well (AUC [ 0.9). However, these differed substantially in prediction and mapped habitat suitability, leading to different ecological understanding and potentially divergent conservation recommendations. Habitat selection was highly scale-dependent and the strongest relationships with environmental variables were at the broadest scales analysed. Modelling approach had a substantial influence in variable importance among top models. Conclusions Our results suggest that optimization of the scale of resource selection is crucial in modelling tiger habitat selection. However, in this analysis,

Introduction
Much of modern ecological analysis is underpinned by ecological niche theory (Hutchinson 1957). At the center of ecological inquiries is the organism, with its distribution, abundance, and survival dictated by its niche, defined as a function of limited resources within an n-dimensional hypervolume of continuous space (Hegel et al. 2010;Blonder 2018). The proliferation of modern analytical tools (Elith et al. 2006;Blonder 2018) and increased computational power to quantify these fundamental and realized niches has advanced the field of landscape ecology and enabled more powerful habitat selection modelling (Hegel et al. 2010). However, such assessments, even for wellstudied species, are rarely straightforward (Mayor et al. 2009).
Gradients of factors along niche hypervolume axes can exhibit complex patterns across a continuum of spatial scales (Cushman 2007;Hegel et al. 2010), with organism responses to its environment is simultaneously influenced by both fine-and broad-scale environmental attributes (Wiens 1976;Johnson et al. 1992;Levin 1992). No single spatial scale is typically sufficient to elucidate organism-habitat relationships. Thus, investigations into habitat selection should be conducted at multiple, ecologically-relevant spatial scales (Wiens 1989;Goodwin and Fahrig 1998;McGarigal et al. 2016). In addition, organism responses to environmental and other factors, may not necessarily be linear. Austin et al. (1990) provide evidence of complex, non-linear organism responses in habitat studies, a result reflected in a number of other studies on species-habitat relationships (Wasserman et al. 2012;Hebblewhite et al. 2014;Timm et al. 2016;Bosco et al. 2019;Macdonald et al. 2019). However, related assumptions pertaining to these relationships have largely been untested.
The importance of incorporating spatial scale into modelling of species-habitat relationships is wellestablished (McGarigal and Cushman 2002;McGarigal et al. 2016). Failing to do so can undermine the performance of habitat selection models and their interpretation (e.g., Wasserman et al. 2012;Mateo-Sanchez et al. 2014;Shirk et al. 2014), potentially leading to errors of inference and application (McGarigal and Cushman 2002). Optimized multiple-scale approaches (see definitions by McGarigal et al. (2016)] have been used to model habitat selection, limiting factors, and threats for a number of species (e.g. Thompson and McGarigal 2002;Toews 2011;Gonthier et al. 2014;Wan et al. 2017), and typically out-perform single-scale models in predictive or explanatory power (Kanagaraj et al. 2011;Wasserman et al. 2012;Mateo-Sanchez et al. 2014;Timm et al. 2016;McGarigal et al. 2016). It is notable, that despite the increased consideration of scale in habitat suitability studies since the publication of the influential McGarigal et al. (2016) review paper, few studies have formally applied scale-optimization to evaluate the effects of landscape components on species habitat suitability at multiple scales.
Conversely, the extent to which the functional shape of a species' response to covariates (i.e. functional form) affects model performance has largely been overlooked. Assuming a single response type for data being modelled may fail to account for different forms of relationship, an assumption that may not be appropriate even for modelling simple ecological associations (Oksanen and Minchin 2002). Indeed, the functional form of species-habitat relationships can vary across habitat characteristics (Redfern et al. 2006) and scales (Wan et al. 2017). Inclusion of variable-specific, non-linear functional relationships may strengthen the relationship between dependent and independent variables, improving representation of covariates in top models that may otherwise be obscured and, potentially, improving explanatory power (Carle 2006).
Generalized linear models (GLMs) have been the most widely-used approach in habitat modelling studies, given their intuitiveness and flexibility . GLMs assume a linear relationship between link function-transformed predictor variables and explanatory variables (Nelder and Wedderburn 1972). However, variables may have unique, non-linear response shapes. The vast majority of habitat relationships studies have assumed linear relationships. In the relatively few cases in which different functional forms have been explicitly investigated, they have typically been restricted to quadratic relationships (Mateo-Sanchez et al. 2014;Mashintonio et al. 2014;Dzialak et al. 2015;Bosco et al. 2019;Macdonald et al. 2019). A few studies have examined other, more complex relationships (Bar-Massada et al. 2011;Fisher et al. 2011;Mateo-Sánchez et al. 2015;DeVoe et al. 2015;Shirk et al. 2018), although this appears not to have widespread application. Importantly, no published study, to our knowledge, has formally compared the performance and interpretation of models with and without functional shapeoptimization.
The effect of scale-and shape-optimization in quantifying species niches and, by extension, specieshabitat relationships merits investigation. Our goal in this paper is to explicitly evaluate the effects of scaleand shape-optimization on modelling habitat selection, using a globally important and understudied tiger population in Thailand as a case study. Tigers are an ideal species in which to investigate this, given large home ranges (Smith 1993;Sunquist 2010), potential influence of broad-scale environmental variation (Krishnamurthy et al. 2016;Reddy et al. 2017), and high conservation priority (Lynam and Nowell 2011;Goodrich et al. 2015). Using long-term camera-trap studies, we evaluate how shape-and scale-optimization affect predictive power and interpretation of models in comparison with non-optimized models. In doing so, we also evaluate how covariates influence habitat selection of tigers and the degree of scale-and shape-dependency.
We test four central hypotheses. First, we expect that habitat selection is highly scale-dependent, with clear patterns of relationship between the scale of analysis and strength of relationship between occurrence and environmental variables. Second, we expect that, for most variables, the relationship with tiger occurrence will be strongest at broad spatial scales, reflecting large home range sizes and high sensitivity to human influences. Third, we expect the multi scaleoptimized model to outperform any of the single-scale optimized models, given evidence from previous studies. Lastly, we expect that shape-optimized models will outperform non-optimized models, given each variable is included with its independent shape of functional relationship with species occurrence.

Study site
The Dong Phayayen-Khao Yai Forest Complex (DPKY) covers 6155 km 2 in eastern Thailand (Fig. 1). It consists of five protected areas-Khao Yai National Park (KYNP), Thap Lan National Park (TLNP), Pang Sida National Park (PSNP), Ta Phraya National Park (TPNP), and Dong Yai Wildlife Sanctuary (DYWS). The complex is relatively mountainous with elevations ranging from 100 m to 1351 m a.s.l. It contains a number of forest types, but primarily consists of mixed evergreen (77%) and mixed dipterocarp/deciduous primary and secondary forest (6%). The complex is biologically diverse, supporting 112 mammalian species, 392 birds, 200 reptiles/amphibians, and 2500 plant species (DNP 2004), and is on the UNESCO World Heritage List in recognition of its outstanding natural and conservation value (UNESCO 2017). Currently, DPKY is almost completely surrounded by villages and human infrastructure.

Camera-trap surveys and data processing
From 2008 to 2017, we conducted camera-trap surveys throughout DPKY to investigate tiger presence in each of its five protected areas (Ash et al. 2020). Survey coverage and intensity increased during this period as funding, camera-traps, and related resources increased, covering diverse habitats across DPKY. Within these areas, cameras were placed to maximize detection of tigers by prioritizing camera placement in areas with previous tiger or prey records and identifying topographic (e.g. ridges, river valleys) or other features (e.g. roads, trails) likely used by tigers (Karanth and Chundawat 2002;Sunarto et al. 2012;Barber-Meyer et al. 2013). To reduce spatial autocorrelation, if more than one camera was present within a 300 m radius, we selected the camera with the highest number of tiger detections or, if no tigers were detected, by greatest survey effort (camera-trap nights [CTN]). A more detailed description of surveys is provided in Ash et al. (2020) and supplementary online materials (Online Resource 1).

Predictor variables
We compiled a set of nine environmental predictor variables based on previous studies of tigers throughout their range (Wibisono et al. 2011;Sunarto et al. 2012;Ngoprasert et al. 2012;Kafley et al. 2016;Thapa and Kelly 2017). These primary variables were further transformed into more biologically informative predictor variables (Table 1) using class and landscape level spatial statistics (McGarigal et al. 2012). Variables were prepared in 30 m resolution to ensure representation of fine-scale variations in the environment.
We used SRTM 1 (Shuttle Radar Topography Mission) arc-second global elevation data (Jarvis et al. 2008), which we further transformed into several variables measuring topographic heterogeneity using the Geomorphometry & Gradient Metrics Toolbox (Evans et al. 2014) in ArcGIS 10.3.1 (Environmental Systems Research Incorporated, ESRI, Redlands, CA, USA, 2011).These included terrain roughness index (degree of elevation difference between adjacent cells), slope position (surface to area ratio) and compound topographic index (CTI). CTI is a measure of flow accumulation (Beven and Kirkby 1979), with low CTI values associated with ridges and mountaintops and high values associated with large, low-lying drainage areas.
Percent forest cover (Hansen et al. 2013) was reclassified into three classes, non-forest (0-20%), open forest (20-40%) and closed forest ([ 40%). Furthermore, the global land cover map (European Space Agency 2015) was classified into eight studyrelevant area classes. Since some forest types (Royal Forestry Department 2000) were similar or underrepresented in our study area, we reclassified eight forest types from an original 12 classes. Streams and rivers were included, but were reclassified to only include perennial rivers/streams (Royal Forestry Department 2000).
To test associations at different spatial scales, we transformed variables into seven multi-scale predictor variables (250 m, 500 m, 1 km, 2 km, 4 km, 8 km, and 16 km). For continuous variables, we generated focal statistics (focal mean and standard deviation) in ArcGIS using a moving circular window with radius equal to each scale. For categorical variables, we calculated landscape metrics in FRAGSTATS (McGarigal et al. 2012) using moving windows for each spatial scale. For class-level statistics, we  Hansen et al. (2013) calculated area-weighted mean radius of gyration (GYR, also known as correlation length), which is a measure of mean potential travel distance within a habitat patch, or habitat extent, and percentage of habitat within the focal landscape (PLAND). At the landscape-level we calculated contrast-weighted edge density (CWED), which quantifies the influence of habitat edges via the degree of contrast between adjacent patches. We also calculated patch density (PD), a measure of the subdivision of habitat within the landscape (i.e. number of patches divided by area). Lastly, we calculated aggregation index (AI), which quantifies the degree of dispersion of habitat patches within the landscape. In total, we prepared 47 variables to investigate how spatial composition and configuration of the landscape at various scales influences tiger detections (Table 1).

Model 1: Single-scale optimized model approach
To provide a baseline for assessing the effects of optimization on model performance, we developed a suite of seven independent models, one for each scale (250 m to 16 km), in which no optimization steps were applied (single-scale optimized [SSO]). This corresponds to what McGarigal et al. (2016) termed pseudo-optimized single-scale (ms4). For these models, the number of variables was reduced through a series of filtering steps. (1) Using the glmer function in  (Burnham and Anderson 2002). These filtering steps resulted in the inclusion of 11 to 19 variables, including survey effort (CTN), among each of the seven scales (Online Resource 2, Table 2). To find the most parsimonious model for each spatial scale we used the 'biostats' package in R (McGarigal 2018), a set of functions which (1) produced models for all combinations of the standardized predictor variables, subsetting the best performing models ranked by DAICc (defined in our models as DAICc \ 2); (2) averaged a final set of models and calculated goodness-of-fit statistics, including proportion deviance explained for each candidate model; and (3) ranked variable importance in the final model with averaged coefficients based on Akaike's model weight (w i ).

Model 2: Multi-scale optimized model approach
For the multi-scale optimized (MSO) model, we selected the scale of best fit for each variable prior to the filtering steps and inclusion in the final model. For each variable, we developed a univariate GLM for each scale and selected the best performing (i.e. most optimal) scale by lowest AICc value. As in other models, we required p \ 0.05 for each variable and selected among correlated variables (Pearson's coefficient |r| [ 0.6) by lowest AICc value. Scale-specific variables were then included in a multivariate GLM using the 'biostats' set of functions (McGarigal 2018) to subset best performing models, calculate goodnessof-fit statistics, and rank variable importance in the final model. This approach is what McGarigal et al. (2016) termed (pseudo-)optimized multiple scales (ms5), as we optimized our model based on a range of scales determined a priori rather than selecting among an unconstrained range of continuous scales.

Model 3: Single-scale and shape-optimized model approach
To investigate the effects of shape-optimization on model performance, we developed seven independent single-scale shape-optimized (SSSO) models, one for each scale (250 m to 16 km), in which variables were included at an optimal functional form. To determine candidate functional forms to investigate, we plotted nonlinear splines using plsmo in the R package Hmisc (Harrell 2018) between variables and binomial tiger detections. Among resulting plots, we identified four potential functional forms (quadratic, logarithmic, exponential, and negative exponential). These functions were used to transform variables at their optimal scales in univariate GLMs, as described above. Unmodified linear forms were also included, resulting in five functional forms tested in this model approach. The best performing shape for each variable was determined by lowest AICc value before conducting, as previously described, filtering based on p-values and collinearity. In the final multivariate GLM for each scale, for quadratic variables, we included both linear and quadratic terms as separate explanatory variables. We were able to include these quadratic terms using 'biostats' (McGarigal 2018), while running all combinations of models, providing an advantage over other similar approaches, such as dredge (Bartoń 2018).

Model 4: Multi-scale-and shape-optimized model approach
In the multi-scale and shape-optimized (MSSO) model, we include both scale-and shape-optimization. As in MSO and SSSO models, we determined the best performing shape and scale for each variable by lowest AICc value before conducting filtering based on pvalues and collinearity. In effect, variables were included in a final multivariate GLM at their optimal scale and functional form. As in the SSSO approach, for quadratic variables, we included both linear and quadratic terms as separate explanatory variables using 'biostats' (McGarigal 2018).

Evaluating model performance and variable importance
We evaluated model performance in three ways. First, we calculated sensitivity (correctly predicted presence/total points present), specificity (correctly predicted absence/total points absent), percent correctly classified (PCC; summary of correctly predicted points at a defined threshold), Kappa (proportion correctly predicted points accounting for probability of chance), and area under the ROC curve (AUC). We measured variable importance via AICc variable importance from the model averaging table and magnitude of the standardized coefficients (Wasserman et al. 2012). We also tested the influence of each variable importance from the model averaging table and magnitude of the standardized coefficients (Wasserman et al. 2012). We also tested the influence of each variable on model performance by calculating the difference in probability of tiger detections when the variable increases from the 10th to 100th percentile, using standardized values and holding all other variables constant at their medians. Variables were back-transformed to original values and plotted to evaluate changes in probability of tiger presence relative to changes in each variable. Lastly, we generated prediction maps to investigate differences in predictions for each modelling approach. These maps depict spatial variation in predicted tiger presence, with each cell containing a predicted presence value ranging from 0 (absent) to 1 (present). To compare predictive maps between the four top models, we used three approaches. First, we calculated the correlation between the pixel values in each of the four maps. Second, we calculated the average absolute difference between the pixel values in the four maps. Third, we classified the four maps into high quality (greater than probability of 0.7) and mid-high quality habitat (greater than probability of 0.5). We then used FRAGSTATS (McGarigal et al. 2012) to calculate percentage of landscape (PLAND), patch density (PD), correlation length (GYR) and largest patch index (LPI) for each of the maps.

Model performance
All models demonstrated high predictive performance, with AUC ranging from a low of 0

Variable scales
Variables were generally represented at broad spatial scales. Almost all variables in the MSSO (six out of seven; Table 3) and MSO (five out of seven; Online Resource 2, Table 4) were selected at the broadest spatial scales (8 and 16 km). In MSO, the smallest scale represented was 1 km, while the smallest scale represented in the MSSO was 4 km. This representation of broad scales is also evident in the similarly strong predictive performance of both single-scale optimized (i.e. SSO and SSSO) 16 km models compared with multi-scale optimized models. It is notable that the worst performing model was the SSO developed at the finest spatial scale (250 m).

Variable sign
The signs of shared variables remained consistent across our four top models. Variables with a positive relationship with tiger occurrence included camera effort (Cam_eff) and percentage of bamboo forest (FT5_PLAND) in all models, standard deviation of compound topographic index (CTI_SD) in the 16 km SSO and MSO models, standard deviation of terrain roughness index (TRI_SD) in the MSSO model, and percentage of evergreen forest (FT1_PLAND) in the

Variable importance
In the 16 km SSO model, AICc variable importance was similar across variables with the number of model subsets in which variables were represented ranging from n = 1 to 2 (Table 4) The 16 km SSSO model exhibited different patterns in AICc variable importance, with the number of model subsets in which variables were represented ranging from n = 1 to 5. Variables with notable standardized coefficients, included percentage of bamboo forest (FT5_PLAND; b = 1.556 ± 0.440), percentage of mosaic cropland (LC2_PLAND; L b = -1.130 ± 0.495; Q b = -0.770 ± 0.516), and camera effort (Cam_eff; b = 1.487 ± 0.311). The strongest effect on probability of tiger presence were evident in camera effort (D0.94) and percentage of bamboo forest (D0.70).
The MSO model had notable variation in AICc variable importance, ranging from n = 1 for percentage shrubland/grassland (LC4_PLAND) to n = 6 for camera effort (Cam_eff) and percentage of bamboo forest (FT5_PLAND). Variable importance, in terms of the magnitude of standardized coefficients, was

Shape of variables
While five functional forms were tested during variable selection for the 16 km SSSO and MSSO models-linear, quadratic, logarithmic, exponential, and negative exponential-only three were represented in fully averaged models. Quadratic relationships (for which linear terms were also included) were evident in percentage of bamboo forest

Model prediction
Correlation between prediction maps (Fig. 3) for the four top models was reasonably high. Correlation of the prediction maps for SSSO-SSO was the highest (0.95), followed by MSSO-MSO (0.92), and MSSO-SSO (0.84; Table 5). Differences between pixel values were generally highest in the central part of the landscape, within TLNP and PSNP, with SSSO and SSO 16 km models having a notably higher predicted presence of tigers over a broader area compared to MSSO and MSO models (Online Resource 3, Fig. 5. Between these models, SSSO had higher predicted values elsewhere in the complex. Predicted tiger presence was higher in a smaller core area of the complex in MSSO and lower predicted presence outside this area compared with other models. The MSO generally had lower predicted tiger presence in more central areas compared to other models though predicted higher probability of presence in certain areas elsewhere in the complex such as in small areas of PSNP and KYNP. Average absolute difference was highest between SSSO and MSO maps (0.17) and lowest between MSSO and MSO maps (0.04; Table 4).
Percentage of high-quality habitat ([ 0.7 predicted tiger presence) in the landscape was highest in the 16 km SSSO (* 12.3%), followed by the MSSO (* 5.5%), the 16 km SSO (* 4.7%), and was notably low in MSO (* 0.7%; Table 6). Extent of mid-high quality habitat ([ 0.5 predicted tiger presence) was also highest in the 16 km SSSO (* 20.8%), followed by the 16 km SSO (* 11.7%), MSSO (* 7.9%), and MSO (* 4.1%). Mid-to high-quality habitat was generally more dispersed in MSO than other models, as indicated by relatively higher PD, and lower LPI and GYR values, characterized by several distinct patches. In the 16 km SSSO/SSO and MSSO, both mid-high and high quality habitat were largely represented by single patches.

Discussion
The focus of this study was to quantitatively evaluate the relative impacts of scale-and functional shapeoptimization on model performance, prediction and interpretation in habitat selection studies, using a globally important tiger population as a case study. Our study confirms the scale-dependence of habitat selection in tigers, and our results suggest it is important for further confirm the importance of habitat selection research to formally evaluate and optimize scale as a component of the modelling process. In contrast, we did not find strong and clear effects of optimizing models to account for different, nonlinear functional shapes of relationship between species occurrence and environmental variables.
Consistent with our first hypothesis, we found that tiger occurrence probability is highly scale dependent. b Fig. 3 Prediction maps of probability of tiger presence, with values ranging from 0 (predicted absent) to 1 (predicted present), generated from the top four fully averaged models-(i) the single-scale optimized (SSO) 16 km model, (ii) the scaleoptimized (MSO) model, the single-scale shape-optimized (SSSO) 16 km model, and (iv) the scale-and shape-optimized (MSSO) model. Map generated with ArcGIS 10.3.1 (Environmental Systems research Incorporated, ESRI, Redlands, CA, USA, 2011)

COR
There were strong differences in the univariate strength of relationship between most variables across scale. In addition, we found large differences in model performance and interpretation, indicating that the ability to predict tiger occurrence is strongly affected by the scale in which variables are measured. This scale-dependence in model interpretation is particularly stark between top models and smaller scale single-scale models (SSO/SSSO). Nine variables in top models were not present in fine-scale single-scale optimized (SSO) models (250 m-4 km) and ten variables represented in these fine-scale SSOs/SSSOs were not present in top models. This is particularly important, as most past tiger habitat modelling, similar to past modelling for other large carnivores (Mateo-Sánchez et al. 2015), have used SSO models with a relatively fine (* 500 m-1 km) focal scale. Our results reinforce that the lens with which the niche of organism is evaluated, specifically its relationships with habitat, has a considerable influence in the interpretation of models. This high degree of scaledependence is evident in other tiger studies (Rostro-García et al. 2016;Krishnamurthy et al. 2016;Reddy et al. 2017), as well as studies on cheetah (Rostro-García et al. 2015), leopard (Pitman et al. 2017;Kittle et al. 2018), puma (Zeller et al. 2014) and clouded leopard Macdonald et al. 2019). Substantial differences in model interpretation and identification of optimal habitat based on the scale of analysis is also consistent with numerous studies of other taxonomic groups (Thompson and McGarigal 2002;Wasserman et al. 2012;Mateo-Sanchez et al. 2014;Wan et al. 2017). Scale-dependence in habitat relationships is a foundational idea in wildlife ecology (Wiens 1989;Levin 1992;Thompson and McGarigal 2002) and critical for analyses which aim to quantify species niches (Hegel et al. 2010). Our results strongly reinforce this and demonstrate the importance of the scale of analysis for predicting tiger presence.
A majority of variables were most strongly associated with tiger occurrence at the broadest spatial scales, which was consistent with our second hypothesis. In the multi-scale optimized (MSO) model, 71% of the included variables were measured at the broadest (16 km) or second broadest (8 km) scale. Similarly, in the multi-scale and shape optimized (MSSO) model, 88% were selected at the broadest or second broadest scale, with more than half at the broadest scale. Notably, there were clear patterns of improvement in AUC in both single-scale model approaches (SSO/SSSO) from fine to broad scales. The 16 km SSO and SSSO models, in which all variables were measured at the broadest scale, were the best performing among these single-scale models. The high consistency among these top models clearly suggests that habitat selection for tigers in our study is largely related to broad-scale habitat factors. These results are aligned with previous studies on tigers (Krishnamurthy et al. 2016;Rostro-García et al. 2016;Reddy et al. 2017), brown bear (Mateo-Sánchez et al. 2015, leopard (Kittle et al. 2018), and clouded leopard , all of which found dominant effects of environmental variables at broad scales. Importantly, however, most habitat relationship studies, even for large bodied and highly mobile animals, typically have utilized single-scale analysis of fine scale (e.g. 500 m-1 km) environmental data . Our analysis suggests that such fine-scale analysis would fail to accurately or fully elucidate species habitat selection patterns, even if model performance is reasonably strong, such as with single, fine-scale models in our study.
Our third hypothesis was that the multi-scale optimized models (MSO and MSSO) would outperform any single-scale optimized (SSO/SSSO) model. This follows from arguments and observations made by McGarigal et al. (2016) who reviewed a large number of habitat modelling papers and found that, in nearly all cases that were formally evaluated, multiscale optimization outperformed single-scale modelling of habitat relationships. The habitat niche of any species consists of a number of dimensions, each composed of environmental variation at particular characteristic scales. Our results are not fully consistent with this third hypothesis. Both the multi-scale optimized and the best single-scale optimized models (16 km SSO/SSSO) performed exceptionally and equivalently well. These results suggest that a single-scale model may sometimes perform as well as a multi-scale optimized model, but only when organism-habitat relationships are dominated by a singlescale response and that the appropriate scale of analysis is evaluated. This has also been observed in other studies (Martin and Fahrig 2012;Krishnamurthy et al. 2016;Timm et al. 2016). In the case of tigers (Krishnamurthy et al. 2016;Reddy et al. 2017), a single, broad-scale model might produce very similar results to a multi-scale optimized model. However, we emphasize that it is impossible to assess the appropriate scale of analyses for a single-scale model a priori. Further, even when a single, large scale may be optimal, broad application can result in potentially important fine-scale niche dynamics being underrepresented (Mateo-Sanchez et al. 2014) leading to oversimplified and potentially misleading interpretations of model results. Therefore, even when a singlescale model provides the strongest prediction, it requires use of scale-optimization to identify and confirm this. Further, when a single-scale model isn't optimal, multi-scale optimization is required to identify the optimal multi-scale model.
Perhaps the most novel aspect of our study is the comparison of the multi-scale and shape-optimized model (MSSO) with multi-scale optimized (MSO) and single-scale optimized (SSO/SSSO) models. We expected that the multi-scale and shape-optimized model and single-scale shape-optimized (SSSO) models would out-perform non-shape-optimized models. This is because, like scale-dependence, species occurrence probability may have a different functional response shape for each variable, with some responses being linear, some unimodal, and some curvilinear (Austin et al. 1990). Results of similar studies suggest that optimizing for functional form may have meaningful influence in the selection of top models (DeVoe et al. 2015;Mateo-Sánchez et al. 2015), potentially allowing for the expression of otherwise obscure environmental effects (Carle 2006), which may improve explanatory power. We hypothesised that optimizing the response shape of variables, like optimizing the scale of response, would measurably improve the predictive ability of the model. Our results do not support this.
While the multi-scale and shape-optimized models were tied as the highest performing model based on the most widely used criterion (AUC), the 16 km singlescale optimized model performed equally well. Furthermore, the multi-scale optimized and the singlescale optimized 16 km models performed slightly better based on other measures. The single-scale optimized (SSO) and single-scale shape-optimized (SSSO) models exhibited similar patterns in model performance, with AUC values generally increasing with scale. While AUC was slightly higher (0.01) in finer-scale shape-optimized models compared to corresponding non-optimized models, this was not the case at broader-scales. This suggests that shapeoptimization may be less important than scaleoptimization.
While model performance was not substantially improved by shape-optimization in our study, the inclusion of variables at their optimal functional form clarified relationships between habitat characteristics and species response. For example, topographic variables (i.e. elevation and terrain roughness) are present exclusively in non-linear forms (i.e. quadratic and logarithmic) in shape-optimized models (SSSO/ MSSO). In effect, the multi-scale and shape-optimized model suggested a more complex relationship between species presence and topographic heterogeneity than in non-optimized models. Expression of these types of relationships are also evident in other species-habitat modelling studies, in particular, the representation of quadratic relationships in final models (Mashintonio et al. 2014;DeVoe et al. 2015;Timm et al. 2016;Macdonald et al. 2018;Macdonald et al. 2019). Notably, inclusion of variables at their optimal functional form influenced representation of variables in final, fully averaged models with differences in variables included between optimized (SSSO) and non-optimized (SSO) models. For example, in 8 and 16 km models, only two variables were included in both SSSO and SSO models out of 12 and 11 total variables in these models, respectively. Further, only two out of 11 total variables were represented in both SSO and MSSO models. This illustrates that models may be sensitive to optimizing for functional form, as with scale, and may lead to considerable differences in model interpretation and conclusions about specieshabitat relationships.

Limitations/caveats
Our fundamental goal in this study was to evaluate the effect of scale-and shape-optimization on the performance, prediction and interpretation of habitat relationships models. We use a single case study species and system (tigers in Thailand) for this evaluation, which limits the generality of our conclusions, however, in exchange we feel our example provides realism in that it reflects actual relationships of a focal species of ecological and conservation importance. We suggest future work employ simulation studies which can control scale and shape effects and, therefore, more reliably and completely evaluate the sensitivity of modelling methods to scale-and shapeoptimization. Second, given our focus was on exploring the novel question of impacts of scale-and shapeoptimization simultaneously, rather than evaluating the best approaches for scale optimization, our example uses GLM modelling because it is the most well understood and widely-used method in habitat modelling, and the dominant approach in multi-scale habitat modelling studies published to date . This approach to scale selection is consistent with many other species-habitat modelling studies (Rostro-García et al. 2016;Kittle et al. 2018;Macdonald et al. 2018). However, with the rapid emergence and adoption of other modelling methods, in particular machine learning approaches (e.g., Evans et al. 2011;Cushman and Wasserman 2018), we suggest future work should focus on comparing performance of scale-and shape-optimization in a wider range of modelling approaches. In particular, tree-based machine learning approaches like random forests perhaps can automatically optimize nonlinear functional shape (Evans et al. 2011), which could be a great advantage, while others, like Maxent, do not (Elith et al. 2011). While there are a myriad of approaches for habitat-selection analysis, we opted to utilize GLMs given their intuitiveness, flexibility, and their wide-spread use in studies investigating scale in similar studies .
We also note that our definition of optimization closely follows the definition proposed by McGarigal et al. (2016) of pseudo-optimization as the scales we evaluated were selected a priori. Importantly, this approach does not optimize scale relationships across a continuous range, and does not optimize them multivariately and simultaneously. A true optimized approach may have allowed for inclusion of more specific scales that best explain variation in tiger occurrence. However, we feel that the scales selected for evaluation were sufficient for understanding relevant scales of relationships between covariates and tiger occurrence and, most importantly, allowed for a more effective and straightforward evaluation of scale-dependence in our study compared to a trueoptimized approach. However, improvement of scaleoptimization approaches is an important current topic in habitat relationships modelling , and we strongly suggest future research to explore the performance of a range of continuous and multivariate simultaneous optimization approaches.

Species management implications
We recommend central DPKY, particularly areas of closed forest cover containing bamboo forest patches within a 16 km window, should be managed as an area of high conservation priority for tigers in this landscape. Our models reinforce the importance of protecting surrounding broad-scale (8 to 16 km) forest within protected areas with low human impact as part of a landscape-scale management strategy. Protection of this core area and facilitating unconstrained movement to other available habitat will be critical to the long-term recovery of this population. Scale-dependence in species-habitat relationships, as demonstrated by this and other studies (Wasserman et al. 2012;Mateo-Sanchez et al. 2014) can have monumental implications for species management. As such, efforts to develop management and recovery strategies should account for scale-dependence of tigers in this landscape. Our results are consistent with other literature highlighting the importance of large protected areas and low human disturbance for tiger populations (Trisurat et al. 2010;Sunarto et al. 2012;Ngoprasert et al. 2012;Hebblewhite et al. 2014;Reddy et al. 2017).
It is pertinent to highlight that our results reflect an expression of the tigers' realized niche within DPKY that may be severely constrained in comparison with its theoretical fundamental niche. The tigers' range itself has suffered considerable declines and its current range likely reflect its refuge in areas less likely to be impacted by human activities, which are the chief driver of range and population declines (Goodrich et al. 2015). In our study area, areas of the landscape with flat topography are typically highly affected by human activities, which may explain their lower selection. Historically, Thailand's lowland forests, now converted by human activity, would have likely represented preferable habitat due to high prey densities (Rabinowitz 1993;Sunquist et al. 1999).

Conclusion
While shape-optimization did not substantially improve performance over other models, it did allow for the expression of potentially important relationships between tigers and covariates that were not apparent in the models assuming linear forms. Importantly, our study clearly reinforces previous studies (McGarigal and Cushman 2002;McGarigal et al. 2016) which highlight the importance of accounting for scale when modelling habitat selection and other ecological relationships, particularly for wide-ranging species such as tigers. We recommend that habitatselection studies on tigers and other species incorporate a robust, ecologically-relevant scale-optimization framework and consider the inclusion and evaluation of shape-optimization in model development.