## Introduction

Spatial distribution of drought plays key role specifically in hydrological research. Drought is a complex natural hazard which can severely affect on a wide range of land, agricultural and socioeconomic infrastructure. Drought can be defined as ‘an extended period of deficient rainfall relative to the multi-year mean for a region’ (Schneider et al., **2011**). Probably, all the climatic zones are suffering from drought, however, its characteristics may vary from zone to zone. Several factors are involved in drought occurrence such as, high wind, low relative humidity, temperature, characteristics and duration of rain, intensity and onset (Wilhite, **1994**). In addition, changed climatic conditions (Dawadi and Ahmad, **2013**) and contaminated water are also causing water scarcity (Aswathanarayana, **2001**). However, demand of water has been increased since past two decades because of enhanced agricultural and industrial activities. Therefore, it is necessary for drought management to define drought conditions for the region under consideration (Wilhite, **2000**).

There are several tools for the assessment of drought, but drought indices are the most commonly used because of its probabilistic structure and universal acceptability around the globe. Drought indices may help in developing effective drought mitigation policies (Kogan, **2000**; Hayes et al., **2004**). Moreover, for the precise and accurate monitoring and quantification of drought, a single climatic variable for drought indices is insufficient (Vicente-Serrano et al., **2010**), as multiple factors are involved in drought phenomena such as, precipitation, humidity, wind speed, temperature, and heat waves (Wilhite, **1994**; Sheffield et al., **2012**; Hao and Singh, **2015**; Vicente-Serrano et al., **2016**). In preliminary stage of the evaluation of drought indices, Palmer et al. (1965) introduced Palmer Drought Severity Index (PDSI) for the characterization and quantification of drought based on precipitation, temperature, latitude and available water capacity. In recent research, it was found that PDSI has the lack of multi-scaling drought study and ability to compare regions (Wang et al., **2017**). Various studies have been conducted to assess the suitability of different drought indices in several applications (Guttman, **1998**; Keyantash and Dracup, **2002**; Musuuza et al., **2016**; Tigkas et al., **2016**; Adnan et al., **2018**). A comprehensive review on the list of available drought indices corresponding with their data requirement is available (Svoboda et al., **2016**).

Parallel to including multiple drought indices in drought monitoring protocols, the role of non-parametric methods have significant impact to explore the information about extreme drought classes (Ghamghami et al., **2017**). In previous research, numerous authors used probability distributions to obtain standardized drought values. Ntale and Gan (**2003**) used Pearson type III (P3) distribution and probability plotting method to improve Standardized Precipitation Index (SPI) values. Blain and Meschiatti (**2015**) provided an alternate procedure to quantify the risk of drought by introducing generalized normal distribution in SPI model. Recently, Stagge et al. (**2015**) proposed multiple distribution for various drought indicators such as SPI and Standardized Precipitation Evapotranspiration Index (SPEI). One of the major advantage of multi-scalar drought indices is to capture drought classification states at various time scales that are essential for both assessing drought in relation to different hydrological systems, and differentiating among different drought types. However, in probabilistic models of each multi-scalar drought indices, uncertainty in accurate and precise estimation of drought classes always exists (Parker, **2014**). Further, the selection of probability distributions in each index is purely subjective in nature. Moreover, it is difficult to include all the time scales in drought monitoring protocol for efficient and accurate drought mitigation policies (Hao and AghaKouchak, **2013**; Hao et al., **2017**).

A non-parametric approach – a graphical method that is based on probability plotting position formula is suggested as a substitute of probability distribution. The fact of using graphical method is to capture the uncertainty of extreme events and hence to reduce errors in accurate and precise estimation of drought indices. Initially, Farahmand and AghaKouchak (**2015**) gave an idea of non-parametric for drought monitoring. They used Gingorten (Gringorten, **1963**) probability position formula as an alternative of Gamma distribution for obtaining standardized drought index. Several other studies used non-parametric approach for drought monitoring (Ghamghami et al., **2017**).

In previous research, various studies used different multivariate techniques to address the multi-scaling issue. Bazrafshan et al. (**2014**) used principle component analysis (PCA) technique to overcome the issue in the selection of appropriate time scale by introducing Multivariate Standardized Precipitation Index (MSPI). Several studies used different time scales to monitor both spatial and temporal behaviour of drought for various climatic regions (Mishra and Desai, **2005**; Patel et al., **2007**; Fiorillo and Guadagno, **2010**). Dogan et al. (**2012**) compared different drought indices using precipitation data with the different time scales and indicated that while comparing drought indices, appropriateness of important time scale has great extent in drought monitoring context. One can quantify a big picture of drought for any region under study without requiring any complex environmental data because of multi-scaling characteristic. Guttman (**1998**) compared different drought indices using precipitation data with different time scales and highlighted that the appropriate selection of time scale has great importance in drought monitoring. Regardless of temporal time scale, classification of affected system due to drought severity for any specified region, spatial scale also contributes significantly (Tsakiris et al., **2013**).

In this study, we aimed to search more relevant and important time scales for Standardized Precipitation Temperature Index (SPTI) (Ali et al., **2017**) at regional level. We proposed Boruta (Kursa and Rudnicki, **2010**) algorithm from the list of machine learning methods for the analysis of 52 meteorological stations of Pakistan. We hypotheses that ranking important time scale using a feature selection approach will lead to overcome the computational difficulties in existing multi-scalar drought indices for future drought monitoring at regional level. Consequently, efficient use of selected drought time scale is possible to inference spatially referenced meteorological stations.

## Methodology

### SPTI – the multi-scaler drought index

There are several procedures to report drought severity using multi-scalar drought index. McKee et al. (**1993**) developed an SPI drought index, which is based on long term precipitation record to quantify the precipitation scarcity. One of the major advantage of SPI index is that it can be used to monitor drought at various time scales. Estimation of quantitative values of SPI can be made by normalizing appropriate probability distributions of the observed monthly cumulative precipitation time series. Consequently, negative and positive SPI values designate less than or greater than median precipitation, respectively.

Parallel to McKee et al. (**1993**) and Vicente-Serrano et al. (**2010**) developed a new multi-scalar drought index: the Standardized Precipitation Evapotranspiration Index (SPEI). In SPEI, water balance model based on the difference between precipitation and potential evapotranspiration (PET) is used with the same mathematical procedure of SPI. One major advantage of SPEI over SPI is to include the effect of evaporation to characterize the regions under study.

Following the same methodology of SPI and SPEI, Ali et al. (**2017**) developed a new multi-scalar drought index; SPTI, to capture drought characterization in both cold and hot climate regions. However, among the list of available drought indices, appropriate selection of particular depends on the availability of data and the nature of climatic zone.

In this paper, we use multi-scalar drought index – SPTI (Ali et al., **2017**) for our case study catchment because of three reasons: (1) the climate of the selected regions has both cold and hot weather, (2) SPTI gives accurate results for regions having low temperature, and (3) SPTI works without any mathematical conflict.

A brief explanation of the stepwise estimation procedure of SPTI is as follows. In the first step, De Marton Aridity Index (DAI) is computed using monthly total precipitation and average monthly temperature for each selected station by the following equation.

Where *P _{i}* is the monthly total precipitation and

*T*is the mean monthly temperature.

_{i}Second step consists on the candidacy of the appropriate probability distribution for *DAI _{i}* series of each station. In this study, 32 most commonly used probability distributions were used to observe the most appropriate probability distribution. . The list of these distributions is available in

*propagate*(Spiess,

**2014**) package in R package.

In third step, fitted distributions are selected for each station’s time series data of *DAI _{i}* that have minimum value of Akaike Information Criteria (AIC) and Bayesian Information Criteria (BIC). For each station, the Cumulative Distribution Function (CDF) of fitted distribution is standardized by using the following method:

In Equation 2, a little modification in the CDF is made to adjust the effect of undefined values in DAI. For example, in case of gamma distribution $q$ is the probability of zero value in the time series data DAI for each station. If $m$ be the number of zeros present in *DAI _{i}* time series data, then $q$ is estimated by $\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$n$}\right.,$ where $n$ is the total number of observation in the

*DAI*time series. Here, instead of the abovestated method, other well-known methods of probability plotting position such as Cunnane (

_{i}**1978**) can be adopted for adjusting probability of undefined values in the CDF.

In last step, the standardized values of SPTI index are estimated by the following approximation procedure which is available in Abramowitz and Stegun (**1965**).

*c*

_{0}= 2.515517,

*c*

_{1}= 0.802853,

*c*

_{2}= 0.010328,

*d*

_{0}= 1.432788,

*d*

_{1}= 0.985269 and

*d*

_{2}= 0.001308 are constant factors (Abramowitz and Stegun,

**1965**).

The resulting standardized values can be classified by using the ranges of SPTI values corresponding to their drought categories as given in Table 1.

Moreover, the estimation of SPTI at various time scales such as 3, 6 and 9 can be made by just taking moving total of *DAI _{i}*series at 3, 6, and 9 months, respectively (Vicente-Serrano,

**2006**).

The selection of appropiriate probability distribution and the calculation procedure of multi-scalar drought index-SPTI at 1, 3, 6, 9, 12, 24 and 48 time scale are estimated by following the guidelines provided in Stagge et al. (**2015**).

### Boruta algorithm and relative importance of time scales

To rank and obtain most important, unbiased, and stable time scales of SPTI drought, a wrapper algorithm boruta (Kursa and Rudnicki, **2010**) is employed using R statistical software. Boruta uses random forest classfier which offers two main adapters: one is permutation importance (raw permutation adapter and normalized permulation adapter) and the other is Gini importance. Gini importance measures the importance of each feature as the sum over the number of splits that include the feature.

Computational strategies of these measures are done using *randomforest* (Liaw and Wiener, **2002**) R package. Additional technical details on the random forest based variable importance measures can be found in Archer and Kimes (**2008**).

Among several machine learning algorithms, the choice of random forest method depends on its two basic reasons: (1) a random forest is the most popular and prevailing ensemble machine learning algorithm and (2) there is no need of manual tuning of its parameter and calculations. We use *Boruta* R package, however several other R packages and software are available (Wright and Ziegler, **2015**).

Random forests include a combination of three predictors, where each tree produces response vectors when it is used as a predictor. In order to estimate error rate and prediction of data, the random forest classifier produces several classifications trees (Liaw and Wiener, **2002**). In Boruta wrapper algorithm with the random classifier, randomness is added to the system and results are obtained from the collective randomized samples by reducing the disingenuous influence of the random variation and correlation. This algorithm calculates weights associated with each feature by iteratively learning the random forest classifier. In our proposed model of Boruta algorithm, spatial reference points of meteorological stations are used as response class, while temporal data on SPTI at different time scales are considered as the predictor variables.

In this research, the steps involved in Boruta algorithm with the random forest as a classifier for the selection of unbiased and stable scale are as follows (Kursa and Rudnicki, **2010**).

- In the first step, for the extension of the range of information system, additional copies of all the variables (time series data on each time scales) are added in the system with a minimum of five shadows.
- In the second step, all the added variables are shuffled to eliminate correlation among each time scale from spatially referenced meteorological stations.
- In the third step,
*Z*-score statistics are obtained by implementing random forest classifier on the entire information system. Where*Z*-score is computed by dividing the average loss by its standard deviation which can further be used as an importance measure. - In the fourth step, a hit is assigned to each time scale by comparing it with those time scales, which have a better score than the maximum value of the
*Z*-score among the shadow attributes. - In the fifth step, a two-sided test is performed for each attribute (variable) with indecisive significance by shadow attributes.
- In the sixth step, permanently removed those original time scales that are significantly unimportant attribute than shadow attribute. In this step, shadow attributes are playing a role as unimportant attributes.
- In the seventh step, attributes which have significant importance over the shadow attributes are then confirmed. In this step, shadow attributes are playing a role as an important attributes.
- In the eight step, all the shadow attributes are removed from the system.

Finally, all the process is repeated for large number times until each time scale in the system receives optimal weights. Hence, the time scales which are statistically significant (i.e. have the higher value of *Z*-score than shadow attributes) are selected as the most important and unbiased time scale. However, those time scales having a small value of *Z*-score than shadow attributes deemed as irrelevant, and are removed from the system.

Before performing Boruta algorithm to rank and search important time scale for socio-geography of Pakistan, we first infer the behaviour of SPTI as a moving average process. The Pearson correlation is used to observe the association among time scales for case study catchment stations.

## Case study

### Data description and study area

We use time series data of monthly total precipitation and monthly mean temperature of $52$ meteorological stations that are located across Pakistan. Pakistan is located in South Asia that is the junction of Central Asia and Middle East, which gives its location great significance. According to the census 2016, its total population is 207,774,520. Where the majority of people belongs to agriculture sectors either directly or indirectly (Rehman et al., **2015**). Its varying regional temperature and precipitation amount in various season indulge the traditional hazard monitoring structure and, require some complex indicators and research frameworks to make exhilarated and effective drought mitigation policies. From past few years, the country suffered severe bad effects which include loss of properties, agricultural lands, prevalence of disease and human deaths, in its historical hazards profile (Spate and Learmonth, **2017**).

In this study, we explore and infer the historical profile of extreme drought events at 52 meteorological stations covering triple of low, moderate and high concentrated climatic variables. Figure 1 shows the selected metrological stations, and further details about selected metrological stations can be found in http://www.pmd.gov.pk/Observatories/index.html.

We use time series data of monthly total precipitation and monthly mean temperature to compute SPTI at vayring time scales. The dataset is collected from Pakistan meteorological department through Karachi data processing centre.

### Calculation of SPTI data

In this study, monthly time series secondary data total precipitation and average temperature from 1955 to 2016 is used to compute SPTI at various time scales by standardizing ${\mathit{DAI}}_{i}$ series. We use Kolmogorov-Smirnov, and Chi-Square, Anderson-Darling tests to find the candidate distribution for each time scale with level of significance $\alpha =.05$ by using Easyfit software (Schittkowski, **2002**). Distributions with minimum value of AIC are then selected for standardization. Further, as different stations have different climatic features, and in each station, the definition of drought with respect to time scale is varied. For example, in moving average process, the definition of meteorological drought at 1- and 3-month time scales, may not lead to monitoring agricultural drought at 6- and 12-month time scale in tropical regions (Dracup et al., **1980**).

Table 2 shows the AIC values of candidate distributions that are fitted on DAI time series at 1, 3, 6, 9, 12, 24 and 48 time scales of Chirat station (Latitude: 33°49.4052′ N, Longitude: 71°53.5752′). At 1-month time scale, 3P Weibull with minimum AIC value –508.52 among all the list of probability functions, have strong candidacy of its standardization. However, generalized extreme values distribution is the best fitted for 3-month time scale, with minimum value of AIC –542.36. Figure 2a and b shows the histograms of those distributions that having minimum values of AIC at all time scale. It is observed that, as the time scale increase, behaviour of data tends to inflate in both sides of tails.

Fig. 2.(a) Theoretical vs. empirical histograms of selected distribution of various time scale of SPTI at Chirat station. (b) Theoretical vs. empirical histograms of selected distribution of various time scale of SPTI at Chirat station.

Figure 2b clearly indicates that the Laplace distribution at 48-month time scale doesn't capture the uncertainty on the significant parts of the data. Moreover, substantial deviation can be seen in extreme right part of the plot. Similarly, Genralized Extreme Values (GEV) at 3-month time scale, Gumbel at 9-month time scale, Johnson at 24-month time scale and Generalized Trapezoidal distribution at both 6- and 12-month time scale are not suitable for drawing accurate inferences. Consequently, the process of selecting the more relevant and appropriate distribution for each time scale may create barrier in fluent drought monitoring system.

Hence, the selection of appropriate time scale for the interpolation of the whole region is the most important task for water management and resources. Moreover, the selection of appropriate drought scale for monitoring drought as compared to the selection of drought indices is more important (McKee et al., **1993**; Vicente-Serrano et al., **2010**). In the next subsection, temporal behaviour of each time scale is assessed using temporal plot and correlation analysis.

### Assesment of temporal behaviour

Figure 3 shows the temporal behaviour of each time scale for Astor, Chirat, Gilgit and Gupis stations. These stations are randomly selected for pictorization of case study computations and analysis. However, similar temporal behaviour of each SPTI time scale was observed for all other stations. On the same line, Fig. 4 presents the correlation and scatters behaviour among each time scale. By careful inspection, it is inferred that as time scale increases, the correlation between indices increases, however, the correlation behaviour decreases from higher to lower time scale. For example, in Chirat station, the correlation between SPTI-48 and SPTI-24 is 0$.73,$*while SPTI-48 and SPTI-1 has very small value 0.15 of the correlation.*

### Selection of time scales

Application of the boruta algorithm with two random forest adapters are applied to confirm and rank the significance of important time scale. The most commonly used time scales SPTI-1, SPTI-3, SPTI-6, SPTI-9, SPTI-12, SPTI-24, and SPTI-48 and spatially reference meteorological stations are used as predictor variables and response class respectively. Results associated with this research show that all time scales are important. However, the SPTI-1 time scale is observed to be the most important by normal permutation and raw permutation adaptor. Figures 5 and 6 show the significant importance of SPTI-1 as compared to other time scales for classification of drought in the selected stations. Hence, SPTI-1 has their own identity and can be used for drought monitoring for spatially referenced geographical studies. Moreover, in both figures the order of relevance, importance of time scales is the same. However, SPTI-9 ranked second by the Gini Index (see Fig. 7). Although the Gini index is a biased method, nevertheless, this method also indicates that all time scales are important.

Therefore, after SPTI-1, it is observed that SPTI-9 in both random forest adapters have a tentative candidate for drought monitoring for this region (see Fig. 7). Higher values of attribute statistics of SPTI-1 confirm its significant importance as compared to other time scales (see Tables 3 and 4). In Table 3, the mean importance of SPTI-1 using random forest raw importance is 0.1558 while median and mode importance is 0.1342 and 0.1398. However, the difference between the minimum and maximum importance is almost negligible.

In Table 4, details about attribute statistics of all scales are described where mean importance of SPTI-1 is significantly higher as compared to other time scales of SPTI. However, Gini index rank 4th to SPTI-1 index (see Fig. 7, Table 5). But, all the time scales are important, and can be used to monitor drought for specific drought type. That is, no single time scale was found unimportant.

## Discussion

This paper presents the solution of issues regarding interpretation and uses of various time scales of multi-scalar drought indices. To overcome multi-scaling issues raised by Bazrafshan et al. (**2014**), this research provides an application of Boruta algorithm for searching and ranking the relevant time scale of the SPTI drought index.

In case study, SPTI at 1, 3, 6, 9, 12, 24 and 48 time scales of 52 maeterological station of Pakistan were computed by standardizing appropriate probability distribution function. To rank and find important time scale of SPTI, Boruta algorithm of machine learning were used to select important time scale. Results of this feature selection approach confirms the significance of all drought time scale. Moreover, experimental results show that SPTI-1, being high ranked and having maximum value of mean importance, determine the exact and joint effect of drought by extracting features from other scales. Two famous and mostly used random forest adapter, suggest the importance of SPTI-1 for spatiotemporal classification of droughts. As different stations have different climatic features, and in each station definition of drought with respect to drought scale are varried. For example, in moving average process, the definition of meteorological drought at 1- and 3-month time scale, may not lead to monitoring agricultural drought at 6- and 12-month time scale in tropical regions (Dracup et al., **1980**). So, the selection of unbiased and stable scale for the interpolation of the whole region is the most important task for water management and resources.

However, for single station monitoring, one can use various time scales depending on the nature of the study. Therefore, instead of the calculation of all time scales for all stations, SPTI-1 is suggested as it is the most important time scale for spatiotemporal analysis. There is no need for the calculation of time series data on all the time scales to infer the regionalized behaviour of drought severity. So, unlike MSPI, one doesn't need to calculate all the other time scales of specified regions. For operational point of view, as each drought time scale has its own importance, so instead of using only SPTI-1, one can use other time scales for a specific purpose. However, for advanced spatial and temporal methodologies, unbiased and stable time scales may contribute in reliable and accurate spatiotemporal interpolations. Hence, for regionalized study and spatiotemporal interpolation of drought, machine learing approach may lead to rank and select the most important time scales of multi-scalar drought indices.

## Conclusions

Current issues in multi-scalar drought indices create several ambiguous problems regarding their computation and interpretation. The purpose of this study is to search more appropriate and unbiased time scale which can be used as a geographical representative, and hence to resolve the issues for accurate and efficient drought monitoring mitigation policies for a particular geography. We present the research framework of how the feature selection wrapper algorithm of Boruta can incorporate for optimal choice of drought time scale in a spatial context. Preliminary, the study used spatio-temporal secondary data on monthly cumulative precipitation and mean temperature. Calculation procedure of multi-scalar drought index-SPTI at 1, 3, 6, 9, 12, 24 and 48 time scale are according the guidelines provided by Stagge et al. (**2015**). Exploratory analysis indicates that as the time scale increase, the correlation between successors and predecessors increase. Therefore, to avoid computational labour and interpretation difficulties of using multiple time scale, appropriate technique of feature selection may be incorporated. In this research, spatially referenced meteorological stations are used as the response variable, whilst all the time scales are playing their role as predictors in Boruta algorithm of supervised machine learning environment. Experimental results show that all the time scales are important in terms of classification. Here, SPTI-1 is the topped ranked and have the ability to capture the spatially referenced variation from other time scales. Similarly, the research may expand for other geographical regions and drought indices. Hence, both the time reduction and multi-scalling interpretation problem can be resolved by using selected drought time scale for inferencing spatially referenced meteorological stations.

## Ethical statement, financial disclosure, and conflict of interest

The manuscript is prepared in accordance with the ethical standards of the responsible committee on human experimentation and with the latest (2008) version of Helsinki Declaration of 1975. The manuscript is prepared by using secondary data and authors have not received any financial support. The authors of manuscript certify that they have no affiliations with or involvement in any organization or entity with any financial interest, or non-financial interest in the subject matter or materials discussed in this manuscript. The authors declare that they have no competing interests.