Estimating Loess Plateau Average Annual Precipitation with Multiple Linear Regression Kriging and Geographically Weighted Regression Kriging

Estimating the spatial distribution of precipitation is an important and challenging task in hydrology, climatology, ecology, and environmental science. In order to generate a highly accurate distribution map of average annual precipitation for the Loess Plateau in China, multiple linear regression Kriging (MLRK) and geographically weighted regression Kriging (GWRK) methods were employed using precipitation data from the period 1980–2010 from 435 meteorological stations. The predictors in regression Kriging were selected by stepwise regression analysis from many auxiliary environmental factors, such as elevation (DEM), normalized difference vegetation index (NDVI), solar radiation, slope, and aspect. All predictor distribution maps had a 500 m spatial resolution. Validation precipitation data from 130 hydrometeorological stations were used to assess the prediction accuracies of the MLRK and GWRK approaches. Results showed that both prediction maps with a 500 m spatial resolution interpolated by MLRK and GWRK had a high accuracy and captured detailed spatial distribution data; however, MLRK produced a lower prediction error and a higher variance explanation than GWRK, although the differences were small, in contrast to conclusions from similar studies.


Introduction
Precipitation, one of the most important climatic factors, is a vital part of the hydrologic cycle, affecting energy transfer and maintaining biosphere functions [1][2][3].It is the focus of hydrology, agriculture, ecology, and environmental science, as well as other related disciplines [4][5][6].Controlling and evaluating the spatial and temporal distribution of precipitation is important in many fields, such as basin and watershed management, soil and water conservation, climate change assessment, agroclimatic or ecological regionalization, ecological environment construction, and prediction and prevention of extreme weather disasters [3,7,8].
The main methods for obtaining precipitation data include ground-based meteorological measurements and spaceborne radar observations [9,10].Spaceborne radar data have low spatial resolution and large uncertainties, which could lead to significant errors in precipitation distribution prediction [11].Ground-based meteorological data are typically used for spatial interpolation of rainfall distribution, due to a longer time series and smaller errors [12].However, as the number and distribution of meteorological stations is limited and inconsistent, it is still a huge challenge to obtain high-accuracy, high-resolution data on the distribution of precipitation [10,13].
Precipitation is a comprehensive reflection of the interactions among the various components of the climate system [2,14].It is influenced by atmospheric circulation, local topography, and land cover, and is closely related to many physical geographic factors, such as altitude, slope, aspect, solar radiation, vegetation type, distance to mountains or seas, lakes, and river systems [15,16].Due to limitations of the interpolation methods algorithm, a lot of studies on this topic have only used elevation (DEM) as an auxiliary variable, combined with co-Kriging or thin plate spline (TPS), in order to predict the distribution of precipitation [6,11,17,18].Many studies have shown that, by considering additional natural geographical features in the interpolation process, rainfall variability can be explained more effectively [5,8,19].Interpolation methods that use a variety of auxiliary variables include multiple linear regression (MLR), geographically weighted regression (GWR), artificial neural networks (ANN), regression Kriging (RK), Bayesian maximum entropy (BME) interpolation, etc. [3,5,8,15,19,20].
Regression Kriging, first coined by [21], is mathematically equivalent to "universal Kriging" (UK) or "Kriging with external drift" (KED) [22].It is a spatial prediction technique, which is combined with a regression forecast of auxiliary variables and Kriging interpolation of the regression residuals.It can combine different regression models to generate many combined methods [20], among which the multiple linear regression Kriging (MLRK) and geographically weighted regression Kriging (GWRK) are the most commonly used.The regression process of MLRK fits with the global trend of the target variables across the study area under stable conditions between spatial variables [5,15].The regression process of GWRK fits the local trend around the predicting points, and can adapt to a non-stationary relationship between the spatial variables, leading to a better explanation of the spatial variation of target variables [5,23].Both of these methods have been widely used in Earth science and environmental science, especially in studies of the spatial distribution of soil properties [23][24][25].There are currently some studies that analyze precipitation interpolation on different temporal and spatial scales using these two methods [5,19,26].However, more research is needed on the use of MLRK and GWRK to evaluate the performances of the two methods and obtain high spatial accuracy precipitation maps.
The Loess Plateau was selected as a study area as there is a shortage of water and severe soil erosion.The region belongs to both semi-humid and semi-arid areas, where vegetation growth is limited by rainfall, and ravines and gullies cover the underlying surface.The relationship between different physical geographical features is significantly non-stationary.Based on the above characteristics, the study area is very suitable for precipitation interpolation experiments.The objectives of this study were as follows: (1) to assess the performance of MLRK and GWRK in the interpolation processes; and (2) to obtain a highly accurate distribution map of average annual precipitation with a 500 m resolution in the Loess Plateau.

Study Area
The Loess Plateau is located in the middle reaches of the Yellow River basin in the north of China, between 33 ˝43 1 N-41 ˝16 1 N and 100 ˝54 1 E-114 ˝33 1 E (Figure 1a).It is surrounded to the east by the North China Plain, to the west by Helan Mountain and the Qinghai-Tibet Plateau, to the north by the Yinshan Mountains, and to the south by the Qinling Mountains, covering an area of 62.29 ˆ10 4 km 2 .The Loess Plateau has a diverse topography.The Luliang Mountains and Taihang Mountains are located in the eastern Loess Plateau and the Liupan Mountains are situated in the western part.The Yellow River flows through the western and northern boundaries, generating the Ningxia Plain and Hetao Plain, where it then crosses the central and southern part of the Loess Plateau, generating the Guanzhong Plain.The middle of the Loess Plateau is covered with highly erodible aeolian silt Water 2016, 8, 266 3 of 20 deposits and has become one of the most severely eroded regions in the world.The northern Loess Plateau contains Mu Us Sandy Land and Kubuqi Desert.A continental monsoon climate is the major climate type of this area.Under its influence, dry and cold winds in winter are followed by frequent and intense rainfall in summer [27].Annual precipitation recorded by meteorological stations ranges from 200 mm to 800 mm and decreases from the southeast to the northwest of the Loess Plateau [28,29].In order to mitigate the prediction error of "edge effects" in interpolation [15], a buffer area with a 100 km bandwidth was considered around the Loess Plateau.
Water 2016, 8, 266 3 of 20 climate is the major climate type of this area.Under its influence, dry and cold winds in winter are followed by frequent and intense rainfall in summer [27].Annual precipitation recorded by meteorological stations ranges from 200 mm to 800 mm and decreases from the southeast to the northwest of the Loess Plateau [28,29].In order to mitigate the prediction error of "edge effects" in interpolation [15], a buffer area with a 100 km bandwidth was considered around the Loess Plateau.

Data Sets and Data Processing
The annual average precipitation distribution in the Loess Plateau was analyzed during the 1981-2010 period.The interpolated precipitation data set was provided by the Climatic Data Center, National Meteorological Information Center, China Meteorological Administration (http://data.cma.cn/).Data came from 277 meteorological stations distributed across the Loess Plateau and 155 stations located in the additional buffer area (Figure 1).
In order to assess the prediction accuracies of the MLRK and GWRK approaches, a validation precipitation data set was selected, including 130 hydrometeorological stations with more than 20 years of data, provided by the National Science and Technology Infrastructure of China, Data Sharing Infrastructure of Earth System Science (http://www.geodata.cn).
In this study, the mean annual maximum normalized difference vegetation index (NDVI), digital elevation model (DEM), solar radiation, slope, and aspect were selected as auxiliary variables for interpolation.NDVI images with 500 m spatial resolution were derived from MODIS MOD13A1 data obtained from NASA EOSDIS Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/).The mean annual maximum NDVI data was the average of the 2000-2014 annual maximum NDVI images (Figure 1b), which were calculated by Maximum Value Composite (MVC).DEM data were collected from the CGIAR-CSI SRTM 90 m database provided by Cold and Arid Regions Sciences Data Center at Lanzhou (http://westdc.westgis.ac.cn) and were combined into 500 m spatial resolution.Solar radiation data were acquired from DEM data using the spatial analyst module of ArcGIS 10.2.Slope and aspect data were also calculated from the DEM data.The unit of slope adopted here was radians.Qualitative aspect data, representing directions, were divided into eight dummy variables where 1 represented data belonging to a certain direction and 0 represented data not belonging to that direction.Directions between 0 ˝-22.5 ˝and 337.5 ˝-360 ˝were assigned to north (N); 22.5 ˝-67.5 ˝to northeast (NE); 67.5 ˝-112.5 ˝to east (E); 112.5 ˝-157.5 ˝to southeast (SE); 157.5 ˝-202.5 ˝to south (S); 202.5 ˝-247.5 ˝to southwest (SW); 247.5 ˝-292.5 ˝to west (W) and 292.5 ˝-337.5 ˝to northeast (NW).In this study, all auxiliary environmental variables were converted into their standardized normal variables, except for the eight dummy aspect variables.

Logit Transformation and Exploratory Data Analysis Methods
Normal distributions of target variables are a requirement for regression analysis and Kriging methods.Figure 2a shows that the Loess Plateau average annual precipitation has an approximately normal distribution.In order to improve this, the logit transformation was applied as follows [22]: where z ++ is logit-transformed precipitation; z + is the precipitation standardized to the 0 to 1 range: z is the precipitation of meteorological station point; and z min and z max are the physical minimum and maximum of precipitation in this region.
precipitation, and the Hot Spot Analysis was used to evaluate the local autocorrelation patterns of precipitation.Furthermore, the logit-transformed precipitation and standardized normal auxiliary variables were used in the Pearson correlation analysis and stepwise linear regression in order to determine the predictors that could apply to the interpolation by regression Kriging.All data analyses and subsequent interpolations were applied using the open-source program of R language [31].

Regression Kriging
Let ( ) i z s represent the target-interpolated variable, where ( ) and n is the number of locations.The estimated values ( ) 0 ẑ s by regression Kriging (RK) can be generally given by Equation (4) [22,32]: One of advantages of using the logit transformation is that the predicted value can be fixed by specifying physical limits (z min ,z max ) based on prior knowledge [22].According to previous studies, average annual rainfall in the study area is generally within 100-1000 mm [29,30].Therefore, z min was set to 100 mm and z max was set to 1000 mm.The logit-transformed precipitation can be reversed to the original scale by using the back transformation: The Moran Index was subsequently used to analyze the global spatial autocorrelation of precipitation, and the Hot Spot Analysis was used to evaluate the local autocorrelation patterns of precipitation.Furthermore, the logit-transformed precipitation and standardized normal auxiliary variables were used in the Pearson correlation analysis and stepwise linear regression in order to determine the predictors that could apply to the interpolation by regression Kriging.All data analyses and subsequent interpolations were applied using the open-source program of R language [31].

Regression Kriging
Let z ps i q represent the target-interpolated variable, where s i pi " 1 ... nq is spatial location and n is the number of locations.The estimated values ẑ ps 0 q by regression Kriging (RK) can be generally given by Equation ( 4) [22,32]: ẑ ps 0 q " m ps 0 q `ê ps 0 q where m ps 0 q is the deterministic part fitted by the regression model, ê ps 0 q is the interpolated residual by ordinary Kriging (OK, a fundamental Kriging method), βk is the estimated deterministic model coefficient ( β0 is the estimated intercept), q k ps 0 q is the value of predictors at the predicted location s 0 , λ i is the ordinary Kriging weight determined by the spatial dependence structure of the residual, and e ps i q is the residual at location s i .
Water 2016, 8, 266 6 of 20 Interpolated residuals by OK can be expressed as: where ˘0 is the vector of OK kriging weights (λ i ), with a constraint ř λ i " 1, and e is the vector of regression residuals.The variance of regression residuals is represented by: where h is the vector of distance and γphq is the semivariance function or variogram.The variogram γphq can be used to solve the vector of Kriging weights λ 0 : where c 0 is the vector of the regression residual covariance at prediction points.

Multiple Linear Regression Kriging
Multiple linear regression Kriging (MLRK) is one of most commonly used regression Kriging methods for predicting the spatial distribution of data.In MLRK, the fitting regression model is multiple linear regression (MLR) with ordinary least squares (OLS) estimates and βk can be expressed as follows: βk " βMLR " ´qT ˆq¯´1 ˆqT ˆz where q is the matrix of predictors and z is the vector of sampled observation.The spatial distribution of logit-transformed precipitation was interpolated by MLRK using six steps: (1) determine the MLR model with predictors selected by stepwise regression; (2) calculate the deterministic part using the MLR model at each prediction point; (3) derive the regression residuals at meteorological stations; (4) filter the optimal variogram model with modelling the covariance structure of regression residuals; (5) interpolate the regression residuals using ordinary Kriging; and (6) add the deterministic part to the interpolated residuals at each prediction point.

Geographically Weighted Regression Kriging
Geographically weighted regression (GWR), which is an extension of multiple linear regression, is also implemented with the ordinary least squares estimates, but takes the spatial locations of data points into consideration.Geographically weighted regression Kriging (GWRK) uses GWR as the fitting model.In GWRK, βk can be expressed as follows: where W is the spatial weighting matrix.The weighting matrix W, which is specified as a continuous and monotonic decreasing function of distance d between the observations, can be calculated by different methods.In this study, the weight of each point was computed by applying the bi-square nearest neighborhood function: where h is the bandwidth for spatially adaptive kernel size.
Water 2016, 8, 266 7 of 20 The optimization of bandwidth was required because a large deviation in regression parameters estimation would be generated if the bandwidth was too large or too small [5,33].The optimal bandwidth can be determined when the cross-validation (CV) scores obtain the minimum value using the cross-validation method [34].The whole procedure of logit-transformed precipitation interpolation by GWRK can be described using seven steps: (1) determine the optimal bandwidth by the cross-validation method; (2) derive the spatial weighting matrix with the bi-square nearest neighborhood function using the predictors selected by stepwise regression; (3) derive the regression residuals at meteorological stations; (4) calculate the deterministic part using the weighting matrix and predictors at each prediction point; (5) filter the optimal variogram model by modelling the covariance structure of regression residuals; (6) interpolate the regression residuals using ordinary Kriging; and (7) add the deterministic part to the interpolated residuals at each prediction point.

Validation Techniques
In order to evaluate the interpolation by MLRK and GWRK, adjusted determination coefficients (adjusted R 2 ) were calculated for the deterministic part by MLR and GWR, and five-fold cross-validation was implemented for interpolated residuals by OK [15].The whole prediction error of back-transformed precipitation corresponding to MLRK and GWRK was evaluated by comparing estimated values with actual observations at validation points.The following indexes were used to verify prediction accuracy: Mean error (ME): Mean relative error (ME r ): Mean absolute error (MAE): ˇˇẑ `sj ˘´z `sj ˘ˇ ( 13) Mean absolute relative error (MAE r ): Root mean square error (RMSE): where l is the number of validation points, z `sj ˘is the precipitation of validation points, and ẑ `sj ˘is the estimation value of the validation points.Adjusted determination coefficients (adjusted R 2 ) were calculated to indicate the amount of variance explained by MLRK and GWRK [5].
where z `sj ˘is the mean of the validation precipitation data.

Spatial Autocorrelation
The closer that two stations are to each other, the stronger the correlation of precipitation. Figure 2b shows that, when the radius of stations was within 25 km, the precipitation correlation of corresponding meteorological stations was as high as 0.909.When the lag was between 100 km to 125 km, the correlation was still above 0.7, indicating a remarkable spatial autocorrelation in this region.In order to quantitatively evaluate the spatial autocorrelation, the Moran Index was calculated, which synthetically measures the autocorrelation for the entire study region.Moran's index I = 0.86 (Z score = 38.50,p < 0.0001), indicated a remarkably significant clustering of precipitation in the Loess Plateau.Furthermore, the local autocorrelation patterns of precipitation can be affected by a non-uniform underlying surface.In order to evaluate the pattern at different meteorological stations, a Hot Spot Analysis was implemented by calculating the Getis-Ord Gi statistics.Figure 3 shows the Z score of local Gi statistics for all considered meteorological stations.It is clear that high Z scores, called hot spots, were clustered in the southern and southeastern Loess Plateau, indicating that precipitation values measured at meteorological stations in these areas were high and similar to each other.On the other hand, stations with low Z scores were clustered in the north and northwest of the Loess Plateau.Z scores of approximately 0, found in the southwestern, central, and eastern parts of this region, indicate an absence of spatial association patterns.
Water 2016, 8, 266 8 of 20 where ( ) j z s is the mean of the validation precipitation data.

Spatial Autocorrelation
The closer that two stations are to each other, the stronger the correlation of precipitation. Figure 2b shows that, when the radius of stations was within 25 km, the precipitation correlation of corresponding meteorological stations was as high as 0.909.When the lag was between 100 km to 125 km, the correlation was still above 0.7, indicating a remarkable spatial autocorrelation in this region.In order to quantitatively evaluate the spatial autocorrelation, the Moran Index was calculated, which synthetically measures the autocorrelation for the entire study region.Moran's index I = 0.86 (Z score = 38.50,p < 0.0001), indicated a remarkably significant clustering of precipitation in the Loess Plateau.Furthermore, the local autocorrelation patterns of precipitation can be affected by a non-uniform underlying surface.In order to evaluate the pattern at different meteorological stations, a Hot Spot Analysis was implemented by calculating the Getis-Ord Gi statistics.Figure 3 shows the Z score of local Gi statistics for all considered meteorological stations.It is clear that high Z scores, called hot spots, were clustered in the southern and southeastern Loess Plateau, indicating that precipitation values measured at meteorological stations in these areas were high and similar to each other.On the other hand, stations with low Z scores were clustered in the north and northwest of the Loess Plateau.Z scores of approximately 0, found in the southwestern, central, and eastern parts of this region, indicate an absence of spatial association patterns.

Exploratory Data Analysis
Correlations between the logit-transformed precipitation (PreT) and the auxiliary variables for each station, including standardized normal variable of DEM (DEM_std), standardized normal variable of NDVI (NDVI_std), standardized normal variable of slope (slope_std), standardized normal variable of solar radiation (rad_std), and eight dummy variables of aspect were tested by

Exploratory Data Analysis
Correlations between the logit-transformed precipitation (PreT) and the auxiliary variables for each station, including standardized normal variable of DEM (DEM_std), standardized normal variable of NDVI (NDVI_std), standardized normal variable of slope (slope_std), standardized normal variable of solar radiation (rad_std), and eight dummy variables of aspect were tested by Pearson correlation analysis (Table 1).Precipitation showed the most significant positive correlation with NDVI, followed by solar radiation, DEM, and slope, where all correlations reached a significance level of 0.01.Among the eight dummy variables of aspect, only the north slope reached a significant correlation level with precipitation. Figure 4a shows that more precipitation (PreT) appears where there is a higher NDVI, and lower precipitation occurs where there is a lower NDVI, thereby indicating that vegetation growth in the region is closely linked to precipitation.With increasing altitude (DEM), PreT gradually decreased and then slightly increased (Figure 4b).The precipitation varied irregularly with a gentle slope, but when the degree of slope increased, precipitation also gradually increased (Figure 4c).The relationship between solar radiation and precipitation is similar to the relationship between DEM and precipitation because of high correlation (r = 0.968, p < 0.01) between DEM and solar radiation (Figure 4d).
Water 2016, 8, 266 9 of 20 correlation level with precipitation. Figure 4a shows that more precipitation (PreT) appears where there is a higher NDVI, and lower precipitation occurs where there is a lower NDVI, thereby indicating that vegetation growth in the region is closely linked to precipitation.With increasing altitude (DEM), PreT gradually decreased and then slightly increased (Figure 4b).The precipitation varied irregularly with a gentle slope, but when the degree of slope increased, precipitation also gradually increased (Figure 4c).The relationship between solar radiation and precipitation is similar to the relationship between DEM and precipitation because of high correlation (r = 0.968, p < 0.01) between DEM and solar radiation (Figure 4d).From Table 1, it becomes evident that there are remarkable correlations between DEM, NDVI, slope, and solar radiation, and correlations between eight dummy variables of aspect also reached a significant level.In order to weaken the impact of co-linearity between the auxiliary variables on regression models, the stepwise linear regression method was used.According to the results of stepwise regression (Table 2), five variables or dummy variables (DEM_std, slope_std, NDVI_std, N, and NW) were selected as predictors and subsequently used in MLRK and GWRK implementation.

Diagnosis and Evaluation of Regression
Based on the results of stepwise regression, the multiple linear regression equation used in this study is expressed as in Equation (17).The significance of regression parameters is given in Table 3. Figure 5a is the scatter plot for the regression-predicted values versus the corresponding residuals at all meteorological stations.The red fitted line near the zero value in this figure indicates that the system correlation between the regression predictions and residuals is weak, which infers a good linear correlation between predictors and precipitation in the regression model.The normal Q-Q diagram (Figure 5b) shows that the regression residuals approximately obey the normal distribution assumption.The scatter randomly distributed around the red fitted line in the scale-location graph (Figure 5c) denotes that the regression residuals were satisfied by homoscedasticity, i.e., residuals corresponding to different magnitudes of predicted values should have the same variance.Figure 5d is the Cook's distance diagram of each station.Cook's distance can artificially reflect the effect of all observation points on the regression model by combining the values of residuals and leverage.The higher the Cook's distance value of a station, the greater the contribution for the regression model.No. 373 and No. 53 meteorological stations, with Cook's distance values greater than 0.06, were regarded as influential observations in this study (Figures 1 and 5d).No. 373 station is located on the edge of Qinling Mountains in the southern Loess Plateau, where the altitude is nearly thousands of metres higher than the surrounding stations.With this sharp rise in mountain altitude, the precipitation was also significantly increased.Northwest Loess Plateau, where the No. 53 station is positioned, belongs to semi-arid areas, where the desert is the main landscape, and rainfall is scarce (100 mm-200 mm).However, due to the influence of Yellow River irrigation on the Hetao Plain, the no.53 station had a higher NDVI than surrounding stations.Thus, the two strong influential points were both caused by local underlying surface changes and could play an important role in accurate simulation of local area precipitation.Geographically weighted regression was also implemented using the selected predictors with the stepwise regression.The optimal bandwidth is 87.8 km in this study because this distance ensured that CV scores reached the minimum value of 45.98 in cross-validation methods.The distribution of local regression coefficients (Figure 6a-e) displayed abundant variation across the study area, especially those of NDVI, north slope aspect (N), and DEM.There were large differences between the mean GWR coefficients and the MLR regression coefficients (Figure 7).All these confirmed the existence of spatial non-stationarity in the study area.The GWR global-adjusted R 2 of 0.93 was much higher than that of MLR (0.44, Table 2), which indicates that the variance explained by GWR was much larger than by MLR.However, the degree of local variance explanations was uneven (Figure 6f) and especially poor in the Taihang and Luliang Mountains.Geographically weighted regression was also implemented using the selected predictors with the stepwise regression.The optimal bandwidth is 87.8 km in this study because this distance ensured that CV scores reached the minimum value of 45.98 in cross-validation methods.The distribution of local regression coefficients (Figure 6a-e) displayed abundant variation across the study area, especially those of NDVI, north slope aspect (N), and DEM.There were large differences between the mean GWR coefficients and the MLR regression coefficients (Figure 7).All these confirmed the existence of spatial non-stationarity in the study area.The GWR global-adjusted R 2 of 0.93 was much higher than that of MLR (0.44, Table 2), which indicates that the variance explained by GWR was much larger than by MLR.However, the degree of local variance explanations was uneven (Figure 6f) and especially poor in the Taihang and Luliang Mountains.Water 2016, 8, 266 14 of 20

Regression Residuals Interpolation
Variograms of multiple linear regression residuals and geographically weighted regression residuals were calculated and then fitted with lag distances.The Ste. Model was found as the most suitable for both residuals (Figure 8).The ratio of nugget and sill in the variogram model (C 0 /(C 0 + C)), called the nugget effect, represents the spatial dependence structure, which can be explained as the proportion of spatial heterogeneity caused by random factors.The higher the ratio, the more variations determined by random factors [35].In general, a ratio of less than 25% indicates that there is a strong spatial dependence structure in the regression residuals; if the ratio was between 25% and 75%, the residuals had a moderately intense spatial dependence structure; and until the ratio reached more than 75%, the spatial dependence structure was very weak, i.e., the regression residuals variability consists of unexplained or random variations.In this study, the nugget effect of MLR residuals was 22.16%, i.e., less than 25% (Table 4), which shows that the spatial variability was predominately caused by structural factors.A nugget effect of GWR residuals of 36.58%signifies that spatial variability caused by random factors in GWR residuals was greater than in MLR residuals.These different degrees of spatial dependence structure in the two regression residuals could be induced by different degrees of trend elimination in residuals of the two regression models.

Regression Residuals Interpolation
Variograms of multiple linear regression residuals and geographically weighted regression residuals were calculated and then fitted with lag distances.The Ste. Model was found as the most suitable for both residuals (Figure 8).The ratio of nugget and sill in the variogram model (C0/(C0 + C)), called the nugget effect, represents the spatial dependence structure, which can be explained as the proportion of spatial heterogeneity caused by random factors.The higher the ratio, the more variations determined by random factors [35].In general, a ratio of less than 25% indicates that there is a strong spatial dependence structure in the regression residuals; if the ratio was between 25% and 75%, the residuals had a moderately intense spatial dependence structure; and until the ratio reached more than 75%, the spatial dependence structure was very weak, i.e., the regression residuals variability consists of unexplained or random variations.In this study, the nugget effect of MLR residuals was 22.16%, i.e., less than 25% (Table 4), which shows that the spatial variability was predominately caused by structural factors.A nugget effect of GWR residuals of 36.58%signifies that spatial variability caused by random factors in GWR residuals was greater than in MLR residuals.These different degrees of spatial dependence structure in the two regression residuals could be induced by different degrees of trend elimination in residuals of the two regression models.In order to quantitatively assess the ordinary Kriging results of two regression residuals, a five-fold cross-validation method was implemented to obtain the prediction residuals, and then the adjusted determination coefficients (adjusted R 2 ) were calculated.The adjusted R 2 of Kriging prediction for the residuals of MLR and GWR corresponded to 0.91 and 0.24, which indicates that the degree of variance explanation of ordinary Kriging was much higher for MLR than for GWR.In order to quantitatively assess the ordinary Kriging results of two regression residuals, a five-fold cross-validation method was implemented to obtain the prediction residuals, and then the adjusted determination coefficients (adjusted R 2 ) were calculated.The adjusted R 2 of Kriging prediction for the residuals of MLR and GWR corresponded to 0.91 and 0.24, which indicates that the degree of variance explanation of ordinary Kriging was much higher for MLR than for GWR.

Validation of MLRK and GWRK
The logit-transformed precipitation predicted by MLRK and GWRK can be calculated by adding the deterministic part by regression to the correspondingly interpolated residuals by Kriging.The logit-transformed precipitations were reversed into the original scale using Equation ( 4), and then the back-transformed precipitation distribution was mapped, as shown in Figure 9.The precipitation distribution, which shows abundant spatial variation, was observably affected by environmental factors such as altitude (DEM) and vegetation.High values of precipitation mainly clustered in the south of the Loess Plateau neighboring the Qinling Mountains where vegetation flourishes.Low values of precipitation appeared in the northwest of the Loess Plateau, which belongs to the sparsely vegetated Ulan Buh Desert and Kubuqi Desert.The back-transformed precipitation of the two prediction methods was statistically analyzed in the range of the Loess Plateau without a 100 km buffer.The range of precipitation predicted by GWRK (100.1 mm-999.5 mm, Figure 9b) is wider than the range predicted by MLRK (136.6 mm-917.3mm, Figure 9a).The precipitation distributions predicted by GWRK showed more spatial variation than those predicted by MLRK.

Validation MLRK and GWRK
The logit-transformed precipitation predicted by MLRK and GWRK can be calculated by adding the deterministic part by regression to the correspondingly interpolated residuals by Kriging.The logit-transformed precipitations were reversed into the original scale using Equation ( 4), and then the back-transformed precipitation distribution was mapped, as shown in Figure 9.The precipitation distribution, which shows abundant spatial variation, was observably affected by environmental factors such as altitude (DEM) and vegetation.High values of precipitation mainly clustered in the south of the Loess Plateau neighboring the Qinling Mountains where vegetation flourishes.Low values of precipitation appeared in the northwest of the Loess Plateau, which belongs to the sparsely vegetated Ulan Buh Desert and Kubuqi Desert.The back-transformed precipitation of the two prediction methods was statistically analyzed in the range of the Loess Plateau without a 100 km buffer.The range of precipitation predicted by GWRK (100.1 mm-999.5 mm, Figure 9b) is wider than the range predicted by MLRK (136.6 mm-917.3mm, Figure 9a).The precipitation distributions predicted by GWRK showed more spatial variation than those predicted by MLRK.Comparing MLRK with GWRK, it is observed that the degrees of variance explanations in two steps of the RK were inconsistent.The adjusted R 2 of the MLR model was substantially less than the geographically weighted regression GWR model; whereas the adjusted R 2 of interpolated Water 2016, 8, 266 regression residuals by Kriging in MLRK was far more than in GWRK.Furthermore, owing to the logit-transformation of precipitation data, the entire prediction error could not be directly evaluated but required complex conversions with the regression models' errors and the ordinary Kriging interpolations [22].Therefore, it was difficult to judge which method ultimately obtained more accurate precipitation predictions.As such, the validation data set was used to calculate the entire errors for two regression Kriging interpolations.The means of back-transformed precipitations (MLRK: 455.3 mm/m 2 ; GWRK: 466.4 mm) were lower than the mean of precipitation validation data (476.2mm).The values of MAE, MAEr, and RMSE of MLRK were better than those of GWRK, but the values of ME and MEr of MLRK were worse than that of GWRK.This implies that MLRK prediction errors on several validation points could be larger than GWRK, but the whole prediction error over all validation points was less than GWRK (Table 5).The degree of variance explanation of MLRK was slightly better than that of GWRK, but in reality, the adjusted determination coefficients (adjusted R 2 ) of the two methods showed little difference.

Discussion
Based on background knowledge of the physical geography of the study area, the geographic factors that are closely related to local precipitation are selected as auxiliary variables.This process can help improve the accuracy of the RK model.The results of this study show that the average annual precipitation in the Loess Plateau gradually decreased from southeast to northwest under the influence of the monsoon and the sea (Figure 9).Affected by complex mountainous terrain, precipitation changed greatly in the Taihang Mountains and Luliang Mountain region in the east of the Loess Plateau.Precipitation in the narrow region windward and to the east of the Taihang Mountains was significantly higher than in the leeward area of the western mountains and the plain area on the east side of the mountains.Due to the Tibetan Plateau, the elevation of the western region of the Loess Plateau gradually increases, and the rainfall increases correspondingly.The Loess Plateau belongs to semi-humid and semi-arid areas, where water availability is an important limiting factor for vegetation growth.Vegetation typically thrives in places where rainfall is abundant.Moreover, if the precipitation changes significantly, the underlying natural geographical factors would also show a marked change.For example, in this study, the reason that station 373 had such a strong influence on observations is that it lies on the Huashan Mountain, where altitude is much higher than the surrounding areas.
By exploring the characteristics of the annual average rainfall data, a suitable interpolation method can be chosen for precipitation predictions.Geostatistical interpolation and non-geostatistical interpolation methods both involve various assumptions.In order to obtain a more accurate prediction, the data need to meet the assumed requirement [22,32].For example, in order to make the distribution of data meet the requirements of the normality assumption, precipitation data were typically adjusted to a logarithmic scale.Of course, it is also possible to use an interpolation method that does not require a normality requirement calculation.
In choosing a suitable interpolation method, it is import to consider whether there is a need for precipitation data to fit trends in the interpolation process, and then further consider whether to fit the global trend or local trends according to the stability of the underlying surface.In this study, MLRK is a global trend-fitting method, which is suitable for a relatively homogeneous underlying surface, and GWRK is a local area trend-fitting method, which is suitable for relatively complex surface types [23,33].Previous studies suggest that accuracy of the local trend-fitting method is generally higher than the overall trend-fitting method [5,19,26]; however, this is not always true.With an increase in the size of the study area and the complexity of the underlying surface, the relationship between rainfall and auxiliary variables will become more unstable.When the underlying surface changes greatly but precipitation is relatively stable, the variability in predicted rainfall may be overestimated in local trend fitting, thus reducing the accuracy of forecasts.This type of situation was recognized in this study.The results showed that the entire GWRK error was slightly larger than the MLRK error.In summary, it is difficult to choose the right interpolation method when starting a precipitation interpolation.The advisable way to obtain better interpolation results is by selecting several available interpolation methods, comparing the error after interpolation, and then selecting the method with the minimum error.
The interpolation result of annual average rainfall is affected by several uncertainty factors: (1) an uneven distribution and limited number of stations leads to the underestimation of trends and rainfall variability in these areas; (2) since the relationship between rainfall and the auxiliary variable is unstable, a more accurate model cannot be precisely established, decreasing the interpolation result accuracy; (3) rainfall station data in most studies are limited to within the study area.A significant prediction error in interpolation appears at the boundary of the study area, known as "edge effects."These effects could be mitigated by taking into account the precipitation data of stations just outside the border [15].In this study, the precipitation data of meteorological stations located at the 100 km peripheral buffer of Loess Plateau was added in order to improve the prediction results of the study area boundary.
Measuring the interpolation error is an important process in interpolation method selection and analysis of interpolation results.In this study, two common methods were used to evaluate the results of interpolation: one using the predictive data set itself with cross-validation; the other using a verification data set to directly calculate ME, RMSE, and other indexes.When using the validation data set, the validation points should cover a broad range of land use types.In this study, the meteorological data of the validation data set is only from hydrological stations at the edge of a river or gully with lush vegetation coverage.This validation data set is therefore lacking in data from other land use types.Therefore, this study can clearly show that MLRK prediction performance was slightly better than GWRK in the area around the river, but in regions with other types of land use, it is not possible to determine which of the two models is better.

Conclusions
One of the main objectives of this study was to generate a highly accurate distribution map of average annual precipitation with a 500 m spatial resolution in the Loess Plateau for the period of 1980-2010.Alternative distribution maps of precipitation were interpolated by MLRK and GWRK methods and all showed a high accuracy.There were large disparities, however, in two regression Kriging processes using different methods: the variance explanation of the GWRK regression model was higher than that of MLRK, but the contrary is true of the Kriging process.The interpolation maps using MLRK and GWRK both captured many details of spatial distribution influenced by predictors.Although the GWRK is based on the spatial non-stationary assumption and the map predicted by GWRK did show greater spatial variation, the final validation analysis revealed that MLRK yielded higher model efficiency than GWRK, with small differences.This is in contrast to other previous precipitation interpolation studies.The conclusions can be summarized as follows: (1) both MLRK and GWRK are able to incorporate multiple auxiliary environmental factors into the modelling process and obtain a highly accurate precipitation distribution map; (2) unlike other studies of precipitation prediction, MLRK is shown to be a better method for precipitation interpolation when the underlying surface is complex.In future, greater effort should be made to consider more physical geographic factors related to or impacted by rainfall as auxiliary environmental variables, which could possess a higher resolution, in order to get a higher accuracy of precipitation distribution maps.More studies should investigate and identify the standards and principles to select and validate the interpolation methods further.

Figure 1 .
Figure 1.Underlying surface features of the Loess Plateau (a); and location of meteorological and hydrologic stations (b).

Figure 1 .
Figure 1.Underlying surface features of the Loess Plateau (a); and location of meteorological and hydrologic stations (b).

Figure 2 .
Figure 2. Structural characteristics of annual average precipitation at meteorological stations.(a) Distribution of annual average precipitation; and (b) correlation of precipitation for different stations at certain lag distance ranges.

Figure 2 .
Figure 2. Structural characteristics of annual average precipitation at meteorological stations.(a) Distribution of annual average precipitation; and (b) correlation of precipitation for different stations at certain lag distance ranges.

Figure 3 .
Figure 3. Mapped Z scores of Gi statistics for meteorological stations.

Figure 3 .
Figure 3. Mapped Z scores of Gi statistics for meteorological stations.

Figure 4 .
Figure 4. Relationship between logit-transformed precipitation (PreT) and main environmental variables, including the standardized normal variable of NDVI (a); the standardized normal variable of DEM (b); the standardized normal variable of Slope (c); and the standardized normal variable of solar radiation (d).

Figure 4 .
Figure 4. Relationship between logit-transformed precipitation (PreT) and main environmental variables, including the standardized normal variable of NDVI (a); the standardized normal variable of DEM (b); the standardized normal variable of Slope (c); and the standardized normal variable of solar radiation (d).
both caused by local underlying surface changes and could play an important role in accurate simulation of local area precipitation.

Figure 5 .
Figure 5. Diagnostic plots for multiple linear regression analysis.(a) is the scatter plot for the regression-predicted values versus the corresponding residuals; (b) is the normal Q-Q diagram; (c) is the scale-location graph; and (d) is the Cook's distance diagram.

Figure 5 .
Figure 5. Diagnostic plots for multiple linear regression analysis.(a) is the scatter plot for the regression-predicted values versus the corresponding residuals; (b) is the normal Q-Q diagram; (c) is the scale-location graph; and (d) is the Cook's distance diagram.

Figure 6 .
Figure 6.Maps of GWR coefficients and local adjusted R 2 .(a) is the distribution of DEM coefficient; (b) is the distribution of slope coefficient; (c) is the distribution of NDVI coefficient; (d) is the distribution of north aspect coefficient; (e) is the distribution of northwest aspect coefficient; and (f) is the distribution of local R 2 .

Figure 7 .
Figure 7.Comparison of regression coefficients between GWR and MLR.The orange points represent the mean of predicator coefficients in GWR, and the blue points represent the regression coefficients of MLR.

Figure 6 . 20 Figure 6 .
Figure 6.Maps of GWR coefficients and local adjusted R 2 .(a) is the distribution of DEM coefficient; (b) is the distribution of slope coefficient; (c) is the distribution of NDVI coefficient; (d) is the distribution of north aspect coefficient; (e) is the distribution of northwest aspect coefficient; and (f) is the distribution of local R 2 .

Figure 7 .
Figure 7.Comparison of regression coefficients between GWR and MLR.The orange points represent the mean of predicator coefficients in GWR, and the blue points represent the regression coefficients of MLR.

Figure 7 .
Figure 7.Comparison of regression coefficients between GWR and MLR.The orange points represent the mean of predicator coefficients in GWR, and the blue points represent the regression coefficients of MLR.

Table 1 .
Pearson correlation matrix between logit-transformed precipitation (PreT) and environmental variables.
** correlation is significant at the 0.01 level (two-tailed); * correlation is significant at the 0.05 level (two-tailed).

Table 2 .
Results of the stepwise linear regression analysis.

Table 3 .
Summary of statistically significant influence of all predictors.

Table 4 .
Variogram models for MLR residuals and GWR residuals.

Table 4 .
Variogram models for MLR residuals and GWR residuals.

Table 5 .
Comparison of MLRK and GWRK performances with validation data.