2.2.1 NR and NRWC Methods

A common method for the imputation of missing data is the Normal Ratio Method (NR), modified by ^{Young (1992)} to the Normal Ratio with Correlation (NRWC). The estimated value V ₀ for the missing data is considered to be a combination of observations with different weights, i.e., V ₀ = (∑ⁿ _i=1 W _i V _i ) / (∑ⁿ _i=1 W _i ), where W _i is the weight of the th nearest meteorological gauge station, W _i is the number of nearby meteorological gauge stations, and V _i is the corresponding observation. Weights for the surrounding stations used in the estimation algorithm were calculated according to Eq. (1):

Where r _i represents the correlation coefficient between the target station and the ith neighboring station, and n _i is the number of points used to calculate the correlation coefficient (^{Xia et al., 1999}; ^{Sattari et al., 2017}; ^{Kanda et al., 2018}).

2.2.2 IDW Method

Another method widely used for estimating missing values in hydrology and climatology is the Inverse Distance Weighting Method (IDW) (^{Radi et al., 2015}; ^{Kamaruzaman et al., 2017}; ^{Sattari et al., 2017}; ^{Kanda et al., 2018}; ^{Morales et al., 2019}). IDW estimation of values, based on an observation, is given by Eq. (2):

where θ _t is the target station estimation, n the number of neighboring stations used in the interpolation, θ _t is the observation at station i, d _it is the distance from the neighboring station i to the target station t, and k is the power parameter, most commonly set to 2, 3 or 4 (^{Teegavarapu and Chandramouli, 2005}; ^{Ford and Quiring, 2014}). IDW has been further modified by several authors, mainly by adjusting the calculation of distance and weighting factors to enhance the outcomes.

2.2.3 CCW Method

^{Teegavarapu and Chandramouli (2005)} proposed improvements to the IDW method by replacing the weighting factor with the correlation coefficient, as follows in Eq. (3):

where r _it is the coefficient of correlation (i.e., the ratio of covariance of two data sets to the product of standard deviations of each data set), derived from all available historical time series data between the data at target station t and their corresponding values recorded at any other neighboring station . ^{Teegavarapu and Chandramouli (2005)} named this new method the Correlation Coefficient Weighting Method (CCW) and showed that it was superior to the traditional IDW method for interpolating missing rainfall values.

2.2.4 MCCW Method

^{Suhaila et al. (2008)} presented several modifications to estimate missing precipitation data from the weighting factors of the NR, IDW, and CCW methods. One of them was the Modification of the Correlation Coefficients Weighting Method (MCCW), in which they changed the weighting function of the CCW method proposed by ^{Teegavarapu and Chandramouli (2005)}, as shown in Eq. (4):

where r _it represents the correlation coefficient of the daily precipitation data between the target station t and the th neighboring station; N is the length of the precipitation time series, and p is a parameter between 2 and 6. Larger values of p assign a greater influence to the values closer to the target station. The most commonly used value for p is 2, according to ^{Suhaila et al. (2008)}.

2.2.5 CIDW Method

Another modification was the Modified Correlation Coefficient with Inverse Distance Weighting Method (CIDW), which is the consequence of combining the IDW and MCCW methods to estimate the missing rain values, as expressed in Eq. (5):

where r ^p _it is the correlation coefficient between the target station t and the ith neighboring station with p between 2 and 6; d _it is the distance between target station t and the ith neighboring station. The minimum distances between the target station and neighboring stations strongly influence the IDW method. Still, the correlation factor could also impact the estimation results, so the proposed CIDW method should improve their imputation (^{Suhaila et al., 2008}).

2.2.6 NRIDW Method

The last of the modifications was the Modified Normal Ratio Method with Inverse Distance Method (NRIDW), which was the consequence of combining the IDW with the NRWC proposed by ^{Young (1992)}. The NRWC method is strongly influenced by positive spatial correlation. At the same time, IDW is affected by the minimum distance between the target station and neighboring stations. Thus, the combination of these weighting factors could improve the outcomes of the estimation of missing values through the weighting factor given by Eq. (6):

The modified methods performed better than the previous versions, according to ^{Suhaila et al. (2008)}.

2.2.7 NRIDC Method

^{Azman et al. (2015)} presented a combination between NR, IDW, and the correlation value, which they called the Inverse Distance Weighting Method of Normal Ratio with Correlation (NRIDC). This method modified the NRIDW proposed by ^{Suhaila et al. (2008)} by including the correlation value. This method keeps the original proposal of combining the normal ratio, the correlation, and the inverse distance in a weighting method, to give more weight to the best estimation to impute missing precipitation data. The weighting of NRIDC is given by Eq. (7):

where r ^p _it is the correlation coefficient between the target station t and the th neighboring station with the best exponent value of p ≥ 4 (^{Azman et al., 2015}); μ _t is the sample mean of the data available at the target station t, μ _i is the sample mean of the data available at the ith neighboring station, and d _it is the distance between the target station t and the ith neighboring station.

2.2.8 HIDW Method

The Altitude Relationship with the Inverse Distance Weighting Method (HIDW) was the consequence of the modification proposed by ^{Golkhatmi et al. (2012)}. They also modified the IDW by inserting elevation parameters into the weighting function, which was optimized using the genetic algorithm. Its weighting function is:

where h _it represents the altitude difference between the target station t and the ith neighboring station, and q and s represent parameters corresponding to distance and elevation, respectively.

2.2.8 GNRIDW and GCIDW Methods

Recently, ^{Morales et al. (2019)} proposed two new generalized weighting methods: Generalization of the Modified Normal Ratio Method with the Inverse Distance Method (GNRIDW) and Generalization of the Modified Correlation Coefficient with the Inverse Distance Weighting Method (GCIDW). GNRIDW constitutes a generalization of the NRWC and IDW methods, in which the weighting factors are as follows:

where, in Eq. (9) and Eq. (10), r _it , d _it and h _it represent the correlation coefficient, distance and altitude difference between the target station t and the ith neighboring station, respectively; n _i is the length of data, q and s represent parameters corresponding to distance and elevation, respectively. GCIDW constitutes a generalization of the CIDW, IDW, and HIDW. It is also distinguished by the inclusion of free parameters and an altitude factor, as in Eq. (10):

^{Morales et al. (2019)} computed the optimal parameters of Eq. (9) and Eq. (10), employing the adaptation strategy of the covariance matrix.

2.2.9 MICE Method

The Multivariate Imputation by Chained Equations (MICE) package for R software (^{R Core Team, 2019}) implements a procedure that allows the imputation of missing records in a given database. It creates multiple imputations (replacement values) for multivariate missing data based on the Fully Conditional Specification (FCS) method, where each incomplete variable is imputed by a separate model. MICE specifies the multiple imputation model based on each variable for a set of conditional densities P (Y _j │X, Y _-j ,R, ϕ _j ), where Y _j is the th column in the data matrix Y, X is a completely observed covariate in the population; Y _-j indicates the complement of Y _j , that is, all columns in Y except Y _j , R is a response indicator of Y, and ϕ _j is a parameter. Starting from an initial imputation, it generates imputations by iterating over the conditional densities (^{van Buuren, 2012}). This algorithm can impute blends of continuous, binary, disordered, categorical, and ordered categorical data. To verify the quality of the imputations, the algorithm generates several diagnostic graphs (^{van Buuren and Groothuis-Oudshoorn, 2011}).

The Predictive Mean Matching (PMM) method is built into the MICE package. This is an attractive way to impute multiple missing data of virtually any pattern, especially non-normally distributed quantitative variables (^{Allison, 2015}).

2.2.10 ReddPrec Method

^{Serrano-Notivoli et al. (2017)} presented the ReddPrec package, developed in R software (^{R Core Team, 2019}), and focused on reconstructing daily precipitation. The methodology incorporated in this package creates daily reference values using all data recorded at the closest stations for each targeted day. To do this, multivariate logistic regression is applied based on the data of the ten nearest neighbors, considering the geographic and topographic variables as covariates. This method optimizes all available information; it does not depend on the length of the precipitation series and preserves the local variability of precipitation distribution.

2.2.11 EM Method

Available in SPSS software (^{IBM Corp, 2017}), the Expectation Maximization (EM) algorithm to fill missing data is an interactive method that estimates unknown parameters of a data model. In applying this method, the missing values are initially calculated using the estimated parameters of the model. The method is based on the reciprocal dependency between the model parameters and the missing values. It consists of two steps: the conditional expectation step and the maximization step. In the first step, the conditional expectations of the missing data are calculated given the observed data and the estimates of the model parameters. In the second step, maximum likelihood estimators of the parameters are found by maximizing the expected likelihood of the first step. These steps are repeated until the iterations converge (^{Schneider, 2001}; ^{Firat et al., 2012}).

2.2.12 RG Method

Available in SPSS software (^{IBM Corp, 2017}), Linear Regression (RG) is a statistical technique that models the relationship between a response variable and one or more input or predictor variables. The result of regression analysis is often the generation of a model that can be used to estimate or predict future values of the response variable. This technique is widely used partly due to the straightforward interpretation of the desired model, and it is commonly used to impute time series of climatological variables (^{Jimenez et al., 2014}). This method requires two steps: one to estimate the relationship between predictors and missing values and another to use a trend equation to fill in the empty data (^{Aieb et al., 2019}).

2.3.1 Performance indicators

Several indicators have been used in the literature to evaluate the performance of missing data imputation methods. These are the Similarity Index (^{Suhaila et al., 2008}), the Variance Ratio (^{Ford and Quiring, 2014}), the Coefficient of Correlation (^{Azman et al., 2015}), the Coefficient of Determination (^{Norazizi and Deni, 2019}), the Mean Absolute Error (MAE) and the Root Mean Square Error (RMSE). The MAE and RMSE have been the most frequently used indicators (^{Ford and Quiring, 2014}). Although both have been used to evaluate method performance for many years, there has yet to be a consensus on the most appropriate metric for evaluating model errors (^{Chai and Draxler, 2014}). RMSE provides a measure of the mean value of the errors of the estimates. Its outcome is in the units as the original observations. However, this indicator should be avoided when there are large measurement errors since these values significantly affect the result when squared. Therefore, RMSE is considered sensitive to outliers, which is why several authors have rejected it (^{Willmott and Matsuura, 2005}; ^{Willmott et al., 2009}). MAE provides an outcome that can be interpreted directly since, like the output of the previous indicator, it has the same units as the original variable (precipitation). Therefore, in this research, we used MAE to compare the performance between different missing data imputation methods since we are using precipitation records without first excluding possible outliers.

The value of the MAE is given by Eq. (11):

where n is the total number of observations, x̑ _i is the estimated value, and x _i is the observed value related to the corresponding meteorological variable in the target station i.

2.3.2 Homogeneity test of data

Homogeneous climatic series can be defined as those influenced only by climate variations (^{Firat et al., 2012}). Several methods have been proposed to analyze the homogeneity of climatological variables.

The Standard Normal Homogeneity Test (SNHT), the Buishand test, and the Pettitt test assume the null hypothesis that annual values Y _i of the test variable Y are independent and identically distributed, as opposed to the alternative hypothesis that a stepwise shift in the mean is present. These tests are able to locate the year where a break is most likely (^{Wijngaard et al., 2003}) to appear. SNHT allows the determination of the inhomogeneous structure at the beginning and/or at the end of the time series (^{Firat et al., 2012}). In contrast, the Buishand range test and the Pettitt test detect a point of change in the observed time series and detect inhomogeneous structures with more sensitivity in the middle of a time series (^{Firat et al., 2012}; ^{Guajardo Panes et al., 2017}). The SNHT and the Buishand range tests assume that the Y _i values are normally distributed, where Y _i (i is the year from 1 to n) is the annual series to be tested, but the Pettitt test does not take into account this consideration (^{Wijngaard et al., 2003}).

SNHT is a likelihood ratio test, and it is performed on a ratio or difference between the candidate station and reference series (^{Peterson et al., 1998}). The comparison of the mean of the first k years of the record with that of the last n - k years is obtained by using the statistic T(k), given by Eq. (12):

where z1-=1k∑i=1kYi-Y-/s, z2-=1n-k∑i=k+1nYi-Y-/s and s is the standard deviation. If a break is located at the year K, then T(k) reaches a maximum near the year k = K. The test statistic T ₀ is defined by T0=max1≤k≤nTk (^{Wijngaard et al., 2003}).

The Buishand range test can be used with variables having any type of distribution (^{Guajardo Panes et al., 2017}). In Eq. (13), the statistic S ^* _k represents the Buishand test, where Y _i (i is the year from 1 to n) is the annual series to be tested, and Y̅ is the mean.

with k = 1,…, n; S ^* ₀ = 0. When the series is homogeneous, the values of S ^* _k will fluctuate around zero, as there are no systematic deviations of the values of Y _i with respect to the mean. If a break is present in the year k = K, then, S ^* _k reaches a maximum (negative shift) or minimum (positive shift) near that year. The significance of the shift can be tested through the adjusted range R of Eq. (14):

where s represents the standard deviation (^{Wijngaard et al., 2003}).

The Pettitt test is a non-parametric rank test (^{Guajardo Panes et al., 2017}) that detects a point of change in the observed time series and is more sensitive for detecting inhomogeneous structures in the middle of a time series (^{Firat et al., 2012}). In Eq. 15, the ranges r ₁ ,…,r _n of each year are used to calculate the statistic X _k of the Pettitt test:

with k = 1,…, n. If a break occurs in the year k = E, then the statistic is maximal or minimal near the year k = E: XE=max1≤k≤nXk (^{Wijngaard et al., 2003}).

^{Wijngaard et al. (2003)} and ^{Guajardo Panes et al. (2017)} classified stations as reliable or useful, moderately reliable or doubtful, and unreliable or suspect according to the number of tests that rejected the null hypothesis. The reliable class includes the stations for which none or, at most, one of the tests rejected the null hypothesis. The moderately reliable class had the stations for which two tests rejected the null hypothesis, and the unreliable class included those for which the three tests rejected the null hypothesis. The results of each test were evaluated for a significance level of 5%.

3.1.1 Semi-arid region: Upper Laja River Basin (CARL)

The MAE values of the 23 stations are in Table III, where the lowest MAE values are identified in bold for each imputation method. The ReddPrec method was optimal at nine stations (~39%); GCIDW was optimal at eight stations (~35%); NR was optimal at three stations (~13%), and NRIDC, HIDW, and EM were optimal at 1 station each (~4%). Overall, the weighting methods (led by the GCIDW method) performed best, with average MAE values ≤ 1.60 mm/day. This was followed closely by the ReddPrec method, with an average MAE of 1.63 mm/day. The rest of the methods had MAE values ≥ 2.28 mm/day. It is worth mentioning that the EM and RG methods imputed negative values, and attempting to force these values to zero should be done with caution, as this can modify the mean and shape of the imputed precipitation distribution. In addition, the MAE values of both methods decreased after forcing negative values to zero, but not enough to change the optimal values presented in Table III.

Figure 2A1 shows the minimum imputed values of selected methods, in which only RG and EM methods imputed negative values. Negative values do not correspond to the physics of precipitation phenomena, resulting in the suggestion by ^{Teegavarapu (2012)} that these methods are not suitable for precipitation.

Fig. 2 LEFT (A1 to A4): Comparison of the model-interpolated values and the observed daily precipitation in the CARL (semi-arid region). The comparison is divided into descriptive values as (A1) minimum (PP_Min), (A2) maximum (PP_Max), (A3) mean (PP_Mean), and (A4) standard deviation (PP_Std) of precipitation. RIGHT (B1 to B3): Comparison of key variables with significant Spearman’s rank correlation in Table IV: (B1) annual mean precipitation and dispersion (Std), (B2) missing data and standard deviation of methods (Std of methods), and (B3) Mean Absolute Error and standard observed precipitation of ReddPrec and GCIDW methods, that were selected as representative methods to simplify the analysis. 

As shown in figure 2A2, there were no relevant differences between maximum imputed and observed values except for station S065, where ReddPrec overpredicted the observed value by 7.9 mm (~10%).

Figure 2A3 shows that the mean precipitation of imputed and observed values had no relevant differences except for the RG method. It is important to note that the mean values of RG and EM methods were evaluated after forcing negative values to zero. Doing so changed the mean annual precipitation from 555 mm/year to 599 mm/year for the RG method and from 553 mm/year to 554 mm/year for the EM method. The RG method generally tended to overpredict the mean observed precipitation in all meteorological stations, with an average difference of ~6.2%.

In addition, in figure 2A4, the standard deviation of the observed and imputed values (negative values of EM and RG methods were forced to zero) among different methods showed appreciable differences at some stations: e.g., S241, S155, S141, S107, S065, S061, S053, S042, and S015. These stations had in common a percentage of missing data greater than 12% (see Table I), except for station S141, which had less than 4% missing data.

Up to this point, it is unclear why certain methods performed better at some climatological stations than others or why the standard deviation of imputed values had larger dispersion between different methods at different stations. We further explored these questions using the GCIDW and ReddPrec methods as representative optimal methods to simplify the analysis.

Spearman’s rank correlation test is a nonparametric approach used to describe the strength and direction of the relationship between two random variables (^{Lyerly, 1952}). In this study, it was applied with key variables and MAE values to identify which variables significantly affected (p<0.05) the performance of GCIDW and ReddPrec methods. Results are provided in Table IV. From this analysis, it is possible to suggest a decrease in methods performance, represented by MAE values, when the dispersion of precipitation data increases (see Fig. 2B3). In addition, the higher the percentage of missing data at the meteorological stations, the greater the dispersion of the standard deviation in the imputed values among different methods (see Fig. 2B2).

Figure 2B1 shows that the increase in mean annual precipitation volume is related to a greater standard deviation of daily precipitation. However, the increase in annual precipitation does not correlate with the rise in elevation.

3.1.2 Humid region: state of Tabasco

The methods EM, GCIDW, MICE, ReddPrec, and RG, were selected to evaluate and compare their performance in a humid region (Tabasco). In this section, auxiliary stations were not considered due to limited availability (see section 2.1).

The MAE values from Table V show that the GCIDW method was optimal in ~59% of the stations (average MAE 6.0 mm/day), EM was optimal in ~24%, and ReddPrec was optimal in ~17%. The RG and MICE methods were not analyzed further because they were not optimal for any station. These numbers confirm a performance reduction of ReddPrec when there are insufficient target stations or auxiliary stations are not considered. However, the EM method (average MAE of 6.5 mm/day) performed better than the ReddPrec method (average MAE of 9.8 mm/day), suggesting that it was not only the inclusion of auxiliary stations that influenced the performance of these methods, but also the precipitation regime.

The EM method predicted negative values (Fig. 3A1) for the semi-arid region. Negative values from the EM and RG methods were forced to zero to evaluate descriptive statistics shown in Figure 3A3 and Figure 3A4. The average MAE value of the EM method decreased to 6.46 mm/day after forcing negative values to zero, but not enough to change the optimal values presented in Table V. Maximum (Fig. 3A2) and mean (Fig. 3A3) precipitation values showed no relevant differences with the observed values, except for standard deviation (Fig. 3A4) at stations T07, T08, T09, T11, T15, T21, T47, T60, T75, T76 and T77, which had a percentage of missing data in the range of 10-24%.

Fig 3 LEFT (A1 to A4): Comparison of the model interpolated values and the observed daily precipitation in Tabasco (humid region). The comparison is divided into descriptive values as (A1) minimum (PP_Min), (A2) maximum (PP_Max), (A3) mean (PP_Mean), and (A4) standard deviation (PP_Std) of precipitation. RIGHT (B1 to B3): Comparison of key variables with significant Spearman’s rank correlation in Table VI: (B1) annual mean precipitation and dispersion (Std), (B2) missing data and standard deviation of methods (Std of methods), and (B3) Mean Absolute Error and standard observed precipitation of ReddPrec and GCIDW methods, which were selected as representative methods to simplify the analysis. 

The ReddPrec method significantly overpredicted the observed maximum daily precipitation at stations T04, T09, T47, T60, and T76 with a difference of 51-220 mm/day. In addition, the observed mean daily precipitation was overpredicted at stations T09 and T47 with values of 1.15 mm/day and 1.44 mm/day (see Fig. 3A2 and Fig. 3A3), as was the observed standard deviation at stations T09 and T47 with an average difference of ~6 mm/day (see Fig. 3A4).

The GCIDW method predicted the observed maximum and mean daily precipitation with no significant differences (see Fig. 3A2 and Fig. 3A3). Nevertheless, there was a considerable difference with the observed standard deviation of precipitation; the standard deviation of precipitation was underpredicted in 45% of stations, with a difference of 0.5-2.5 mm/day.

The effect of key variables on the performance of the methods is summarized in Table VI, where again, the percentage of missing data has a significant correlation with the dispersion of GCIDW, EM, and ReddPrec standard deviation (Fig. 3B2). Elevation did not have a significant correlation with the MAE values; the explanation is probably due to the low height of the stations, which varies between 2 m and 83 m with an average of 23 m (see Table II).

Annual precipitation volume was significantly positively correlated with the standard deviation of daily precipitation (Fig. 3B1). However, it is impossible to suggest a decrease in performance, represented by MAE values, with increasing precipitation dispersion and missing data percentage (see Fig. 3B3).

3.1.3 Precipitation regime influence over MAE values in the semi-arid and humid regions

Given the results analyzed in the semi-arid and humid regions, where key variables with significant correlation were analyzed separately by climate region, it is also relevant to analyze both climate regions together to show the evolution of MAE values for GCIDW, ReddPrec, EM, MICE, and RG methods in function of the mean annual precipitation of each meteorological station (see Fig. 4).

Fig 4 Comparison of the MAE values for selected methods and the annual mean precipitation of meteorological stations in the semi-arid region (CARL) and the humid region (Tabasco). 

In Figure 4, it is possible to observe that the performance of the GCIDW and EM methods gradually decreased as the precipitation regime increased, up to a threshold of ~2000 mm/year, at which the performance of these methods stabilizes with a slight tendency to improve. On the other hand, the RG and MICE methods showed a clear tendency to decrease in performance as the annual mean precipitation increases. Finally, there was a significant correlation between the MAE value of the ReddPrec method and the annual mean precipitation (see Table IV and Table VI); however, this linear correlation is less evident in the humid region than in the semi-arid region, especially in the range between 1500-2500 mm/year, where MAE values are higher than 14 mm/day at stations T09, T15, T34 and T37 which visually do not follow a linear correlation (see Fig. 4).

3.2.1 Semi-arid region: Upper Laja River Basin (CARL)

Table VII summarizes the p-values and the years of four stations (17.4%) in which the homogeneity tests detected changes in precipitation, including two stations (S061 and S141) where three homogeneity methods coincided in detecting the year 2001 as the year of change. The year 2001 was also the most frequent in terms of changes in homogeneity. This change in the homogeneity could be related to the modernization of the equipment and instrumentation of the meteorological stations from 2001 to 2006 (^{CONAGUA, 2012}). 78.3% of stations were homogeneous, without changes in precipitation. The homogeneity test was applied individually for 14 imputation methods. The results showed the same number of reliable stations for all 14 but with some differences concerning the year of change in the unreliable or moderately reliable stations.

3.2.2 Humid region: state of Tabasco

Table VIII shows p-values and the years of three stations (10.3%) in which the homogeneity tests detected changes in precipitation. Twenty-six stations (89.7%) were classified as reliable, with the precipitation data corresponding to homogeneous conditions.