Controlled Sampling Approach in Improving Multiple Imputation for Missing Seasonal Rainfall Data

Missingness in rainfall data is one of the well-known and challenging issues faced by meteorologists and researchers from all over the world. The problem would affect the quality of the data which is very important in representing the actual meteorological characteristics of a particular location. Therefore, the missing data should be properly treated in order to provide good quality dataset for the public domain. In furtherance of ensuring the accuracy of imputed missing data, the original structure of the rainfall data series must be specifically preserved when the data are having seasonal patterns. Most of the environmental datasets are generally characterized by outliers and seasonal patterns. These characteristics have certainly affected the performance of missing data imputation methods. The problem of missing data can be treated, but a specific structured approach must be employed when involving dataset that contains outliers and seasonal patterns. This study has highlighted and discussed the structured and comprehensive procedures on how to tackle the problem of missing data by emphasizing on controlled sampling approach for their implementation. The missing values were estimated by using multiple imputation based on block bootstrap approach associated with normal ratio methods compared to the conventional sampling (i.e. general bootstrap approach). The analysis and experimentation are illustrated using several datasets obtained for several locations in Peninsular Malaysia. The block bootstrap approach has revealed its advantage of preserving time series structure in its process and successfully improved the estimates of missing rainfall data imputation.


Introduction
Quality environmental data are highly necessary for performing efficient and effective hydrological and meteorological analyses. Completeness in the data is one of the most important factors to be considered in ensuring the quality of any data. However, the occurrence of missing values in environmental data, especially rainfall data is inevitable (Di Piazza et al., 2011;Miró et al., 2017) and has been a long-standing problem that needed to be addressed (Ramos-Calzado et al., 2008). There are various reasons that contribute to the missingness of rainfall data. Factors such as meteorological extremes, observational recording errors, relocation of stations, malfunctioning of instruments, and human errors during the measurement of rainfall data may affect the consistency in the rainfall records.
Missing data is considered one of the barriers in the production of accurate results in rainfall analysis. Missing data in the process of analysis may result in biasness and produce misleading results (Ramos-Calzado et al., 2008;Kalteh and Hjorth, 2009). As a result, inaccurate information will be passed on to the hydrological management and development activities such as water resource management, irrigation scheduling, and flood forecasting and prediction. Therefore, this shortcoming has motivated researchers to identify suitable imputation approaches in treating missing data.
A common practice when dealing with the missing values is to simply discard them from the rainfall records in the data preprocessing stage (Acuña and Rodriguez, 2004). However, it is not recommended because significant information could be lost when incomplete data is neglected. This will subsequently cause a loss of efficiency and produce severely biased estimates (Donders et al., 2006;Aydilek and Arslan, 2013). Therefore, to ensure that an effective rainfall analysis is performed, it is essential to impute the incomplete rainfall dataset by using appropriate methods.
It is important to understand the reasons and mechanisms behind the missing data for the observed dataset before introducing the missing values in the dataset. Generally, the missing data is introduced based on three fundamental mechanisms, i.e. i. missing completely at random (MCAR), ii. missing at random (MAR), and iii. missing not at random (MNAR) (Twumasi-Ankrah, Odoi, Pels, & Gyamfi, 2019). The missing mechanism is categorized as MAR if the probability of having missing data in a variable depends on observed data, however, it is known as MNAR if the missing data are related to unobserved data. On the other hand, the mechanism is called MCAR if the missing data is independent of both unobserved and observed data. These mechanisms are essentially related to the underlying reason why the data is missing and explain the relationships between measured variables and the probability of missing data (Jahan et al., 2018).
Meanwhile, it is necessary to identify appropriate treatment approaches for missing values based on the mechanism (or pattern) of missing data (Hanaish et al., 2013;Jahan et al., 2018). According to Kalteh and Hjorth (2009) and Yozgatligil et al. (2013), MCAR and MAR are called ignorable response mechanisms because the reasons for the missing data can be ignored during the analysis. Model-based methods require this assumption, and it is reasonable to assume that missing hydrological data are MCAR or MAR (Kalteh and Hjorth, 2009). Meanwhile, the NMAR is a nonignorable pattern of missingness. Ingsrisawang and Potawee (2012) and Jahan et al. (2018) assumed that the existence of missing rainfall values in their studies were under the MAR mechanism. On the other hand, Hanaish et al. (2013) claimed that the missing values in Malaysian rainfall data is assumed to be in MCAR category, which means that the cause of missing data is either not related to the values that are missing or is not dependent on the observed data, totally random missingness (Fielding et al., 2009). In a similar study by Farhangfar et al. (2004), the missing values under MCAR mechanism was introduced in their study. From the reviewed literature, it can be concluded that most of the missing values in hydrological data are categorised under MAR mechanism because the presence of missing data only depends on the observed data rather than the other variables (e.g. the latitude, the longitude, or the elevation of the rainfall stations) (Chuan et al., 2019). Since several missing data estimation methods are to be used in this study, including multiple imputation approach, the missing data mechanism can be assumed to be missing at random (MAR).
Various approaches have been suggested and developed by previous studies for estimating the missing rainfall values. However, multiple imputation (MI) is the most popular treatment among the other approaches used in the last few decades due to its high capability in handling the missing data. Several studied have found this approach to be a powerful tool in most of the situations (Yendra and Jemain, 2013;Yozgatligil et al., 2013;De Carvalho et al., 2017;Ekeu-Wei et al., 2018;Chen et al., 2019). Multiple imputation approach was introduced by Rubin (1987) to overcome the limitations of single imputation (SI) approach. This approach provides a more accurate and robust estimation results by considering the uncertainty of missing values and variability of the imputed values in the imputation process (Lo Presti et al., 2010). Multiple imputation approach treats the missing data in a way to produce a valid statistical inference instead of estimating the missing values as close as possible to the observed ones (Enders, 2010). Enders (2010) organized the MI approach into three basic phases, i.e. the imputation phase, analysis phase, and pooling phase. Generally, imputation is the most important phase in the MI approach. It involved the imputation of the missing values as in the process of SI approach. The process of imputation is implemented multiple times (m times); in many cases, three to five times (Rubin, 1987) on different datasets. The imputed values are different for each dataset due to the inherent randomness in the algorithm itself. The complete imputed datasets are then analysed by standard methods which result in m sets of parameter estimates. Finally, these estimates are pooled to produce a single point estimate. Yozgatligil et al. (2013) and Jahan et al. (2018) adopted the Monte Carlo Markov Chain based on the expectation-maximization (EM-MCMC) method in the application of MI approach in their studies. The method considers missingness as a proportional information of the sample to estimate the parameter of interest through conditional expectations (Yozgatligil et al., 2013). The MCMC estimates the imputation and parameters iteratively with the initial estimates provided by the EM algorithm. The EM-MCMC method has been found to produce robust estimation results, nevertheless, it was unfortunately proven by Yozgatligil et al. (2013) as inefficient in the application of MI approach. It takes longer run times with more intensive computations compared to simple conventional methods such as the arithmetic average and normal ratio.
In the last few decades, MI approach was frequently applied by several researchers (Khalifeloo et al. (2015), De Carvalho et al. (2017, Miró et al. (2017), Sattari et al. (2017), Jakhar et al. (2018), and Milo et al. (2019)) in the imputation of missing rainfall data. Due to its efficiency, there are various MI approach packages that have been introduced for missing data problems. R is one of the programming languages that offer MI packages. The packages comprise of "MICE", "Amelia II package", "missForest", "Hmisc", and "mi" however, the Amelia II package (Amelia) is the one that is mostly used in imputing the missing rainfall data.
Amelia is a bootstrapping-based algorithm with the adoption of the expectation maximization (EM) method that is applied in the estimation process (Honaker et al., 2018). Bootstrap is a resampling technique that is used to estimate statistics on a population by sampling a dataset with replacement (Bland and Altman, 2015). The bootstrap creates multiple datasets from the observed dataset without the need -to make any assumptions. The created sample is a set of randomly chosen observations that has the same size and equally representative of the original observed dataset. This was the main reason for considering bootstrap approach in multiple imputation. The bootstrap approach allows for the occurrence of uncertainties in producing the estimated values. The general bootstrap is the conventional approach that performs generalized random sampling. General bootstrap was considered in Amelia to simulate the estimation uncertainties in implementing multiple imputation. Cardenas et al. (2009) revealed that the weakness of the general bootstrap is that, it does not preserve the time series original structure in its process, and this makes the conventional sampling approach unsuitable to be used for data that has seasonal patterns, such as rainfall data.
According to Yunus et al. (2011), rainfall distribution in Peninsular Malaysia is commonly governed by two rainy seasons, i.e. southwest monsoon which usually occurs in mid of May and ends in August and northeast monsoon which initiates in early November and ends in February. There are transitional periods between the monsoon seasons which usually occur in April and October. The northeast monsoon season brings heavy rain to the east coast areas (involving the southern and eastern parts) of Peninsular Malaysia (Jamaludin et al., 2010) while the west coast region experienced heavy rainfall during the southwest monsoon season ). Information about the seasonal element is very important to be considered in dealing with the rainfall dataset (Deni et al., 2009). The original structure of the rainfall time series must be preserved so that all the significant information can be used to ensure accurate results of the analysis being conducted (De Carvalho et al., 2016).
Therefore, block bootstrap with controlled sampling technique was introduced in the MI approach to improving its performance, especially when dealing with rainfall time series data. This study is aimed to introduce the block bootstrap for data sampling purposes in the MI approach. This is an effort to improve the performance of the currently existing approach which uses generalized random sampling (general bootstrap) in order to produce accurate imputed values in providing a good quality dataset to be used for the public domain.

Data description and study area
The daily rainfall data was collected from the Malaysian Drainage and Irrigation Department and the Malaysian Meteorology Department. It comprises of data from 13 rainfall stations throughout the West and Southern regions of Peninsular Malaysia for the 40-year period between 1975 to 2014. Most of the rainfall stations recorded data using automatic tipping bucket rain gauge, which has a sensitivity of 0.5 mm per tip. However, for some stations, the data collection is still executed using manual methods of measurement. The amount of rainfall for a particular day is the amount collected over the 24-hour period starting at 8:00 a.m.
Lalang Sg Lui (West region) and Johor Bahru (Southern region) were chosen as the target stations due to the following reasons: (i) provide supply to water resource (i) existence of main economic and administrative centers and (ii) serves as the central industrial and commercial hub. The neighbouring stations are chosen by considering the radius distance to the target station as followed by Jamaludin et al. (2008) and Jahan et al. (2018). The stations located within the radius of 100 km from the target station were selected as the neighbouring stations. This is the optimal distance as it considers a suitable number of stations and gives reasonable estimation results (Jamaludin et al., 2008). A greater radius could possibly slow down the computation time as the number of stations increases while a smaller radius will result in no neighbouring station available within the range. Therefore, 100 km was chosen as the best distance for selecting the neighbouring stations. This is due to the same study area (Peninsular Malaysia) considered in both studies which are involved with almost the same distribution of rainfall stations. Three neighbours were assigned to Johor Bahru while eight neighbours were chosen for Lalang Sg Lui station. All the considered stations with their respective geographical coordinates and spatial information are listed in Table 1. Their locations in the geographical map of Peninsular Malaysia are presented in Fig. 1.

Table 1
The target stations (in bolded) and their neighbouring stations within the respective region with geographical coordinates and distances.  2 illustrates the stages involved in missing data imputation analysis. The multiple imputation procedure starts with the process of resampling for the observed dataset. Two sampling approaches are considered, i.e. general bootstrap which represents the conventional approach and block bootstrap denotes the controlled sampling approach. Then, the missing values are estimated using several normal ratio methods such as Old normal ratio (ONR), Modified normal ratio based on the trimmed mean (NRTR), Modified normal ratio based on the median (NRMED), and Modified normal ratio based on the geometric median (NRGMED). Finally, the performance of the estimation methods for both of the sampling approaches are compared based on several criteria, namely, MSE, RMSE, MAE, and SIndex. Their performances are evaluated based on various percentages of missing data and the outliers that will be created in the dataset. This way, the accuracy and consistency of the results is spontaneously assessed. The detailed explanation for each stage will be discussed in the next section. The experimental design executed in this study for the imputation of missing values in the daily rainfall dataset is discussed below. The procedure for the experimentation is conducted based on several stages, i.e. data handling, sampling approach, missing values imputation, and performance evaluation.

Data handling
There are two different types of procedures employed in the study. The first procedure utilizes the observed dataset whereas the second procedure uses the same dataset with few numbers of created outliers. The datasets used in the first and second procedures are named Dataset 1 and Dataset 2, respectively. Fig. 3 presents the data handling procedure implemented in this study.
In order to assess the robustness of the imputation approaches in this study, outliers were created in the observed dataset of the target station, which is named Dataset 2. The outliers were created based on three different levels, in terms of outliers' percent (namely 5%, 10%, and 15%). The procedure for the creation of outliers to be used in the rainfall dataset is explained below: i. Set the number of outliers to be included in the dataset for the experimentation (i.e. 5%, 10%, and 15% of the dataset). ii. Identify the existing outliers in the observed dataset of the target station using the boxplot method as applied by Laurikkala et al. (2000). In the boxplot, the boundaries are set to determine the outlying observations. The observations that exceed the boundaries are considered as outliers. iii. Determine the percentage of existing outliers in the observed dataset. iv. Create outliers that have been set in (i) by considering the percentage of available outliers in (iii).
The artificial outliers are generated from the uniform distribution as suggested by (Tax and Duin, 2001). v. The generated outliers are randomly replaced with the observed values in the dataset. Since a complete observed rainfall dataset was used for this study, the missing data were artificially created in both datasets (Dataset 1 and Dataset 2) in order to evaluate the performance of the imputation approaches. Three different levels of missingness percentage (%) were considered to assess the consistency of the estimation results. The missing values are imputed by the considered estimation methods using these following procedures: i. Identify a target station that has problems of missing values, however, for the purpose of analysis, the station with a complete dataset is used. ii. Identify neighbouring stations based on a 100 km radius distance from the target station. iii. Resample rainfall dataset of the identified target station (i) to produce bootstrapped dataset. iv. Repeat step (iii) five times to produce five bootstrapped datasets. v. Generate artificial missing values in the bootstrapped datasets in (iii) for the three different levels of missingness percentage (5%, 10%, and 15%) randomly.
(Supposedly, if 10% of the data that is randomly chosen for testing methods are considered missing, the rest 90% of the data are used to calculate the estimation values for missing data).
vi. Using the dataset of neighbouring stations obtained in (ii) to estimate the missing values in datasets in (iv) using the proposed estimation methods. vii. The missing data is then imputed with the estimated values obtained in (vi) to produce the complete rainfall datasets. viii. Do error analysis on each of the complete imputed dataset.
ix. Average the results of the analysis to produce a single result.
x. The results in (ix) are then compared between the estimation methods to evaluate their performance.

Sampling approach
In this study, two MI approaches were considered, i.e. MI based on general bootstrap (conventional sampling) and MI based on block bootstrap (controlled sampling). The former MI approach is named MIboot while the latter is named MI-blockboot. Table 2 represents the comparison of both MI approaches.
The MI-blockboot is the enhancement of the conventional sampling MI approach. The block bootstrap was introduced to overcome the limitation of general bootstrap incorporated in the MI-boot wherein the original structure of the rainfall series is not preserved. Based on the comparison, the MI-blockboot shows its advantages over the MI-boot. It's ability to preserve the original structure of the time series dataset promises more accurate estimation results. Block bootstrap (also known as moving block bootstrapping (MBB)) was exclusively introduced for time series data (Inoue & Shintani, 2006). This is due to the characteristics of time series dataset wherein there is a need to preserve its original structure. Therefore, in order to make sure that the original structure of the studied dataset is not affected, block bootstrap was introduced in the current study. The adoption of block bootstrap in the MI was intended to improve the results of imputation, especially when dealing with time series data. The block bootstrap was pioneered by Carlstein (1986) and further developed by Kunsch (1989). The idea featured is to resample blocks (intervals of a time series) with replacement. The blocks are usually disjointed and cover the entire time series, but they can be overlapping (Chernick and LaBudde, 2011). Block bootstrap divides the rainfall time series into several blocks and samples all of the blocks before concatenating them. As a result, the dependency structure of the time series was preserved within each block (Li, 2006). Block bootstrap was adopted by Burhanuddin et al. (2017b) in their MI that was specifically proposed for rainfall time series. The bootstrapping approach was applied to normal ratio methods to impute Data are divided into several blocks before being resampled Data are directly resampled missing values in the Malaysian daily rainfall dataset. They have proved the capability of the proposed MI in providing accurate estimation values, especially when dealing with the dataset that contains outliers. The procedure was a bit different from other studies in determining the size and number of blocks. Rainfall time series were divided into several blocks according to the characteristics of the rainfall data, such as the seasonal pattern and the rainfall amounts before going through the actual bootstrap process.
The bootstrap was applied separately on each block of the data. The detailed procedure of block bootstrap that has been implemented in the current study is explained below (Burhanuddin et al., 2017b).
i. Divide rainfall time series into several blocks.
Rainfall time series distribution is investigated based on the seasonal pattern. Monsoon seasons in Malaysia are reviewed to understand the pattern of rainfall amounts throughout a certain period of a year. Generally, the monsoon season is divided into three categories, i.e. southwest monsoon, northeast monsoon, and inter-monsoon. The southwest monsoon is the source of a high total amount of rainfall in the northwest, southwest, and west regions of Peninsular Malaysia which usually occurs from the middle of May and ends in August. On the other hand, the northeast monsoon which usually begins in early November and ends in February only caters to the eastern region of Peninsular Malaysia with a high total amount of rainfall. Inter monsoon period occurring in between the two monsoons, which are during the March/ April and September/ October is also usually associated with a heavy rainfall. Therefore, from the reviews, the rainfall time series is divided into four blocks, i. ii. Resample each block of rainfall time series obtained in (i) for 1000 times (the size of the bootstrapped sample is the same size as the original data). Then merge the blocks of the bootstrapped sample to produce daily rainfall time series for a year. iv. Repeat step (i) to (iii) for 5 times to produce 5 new samples The new samples obtained from the bootstrapping process were used in the next imputation procedure discussed in detail in the following section.

The imputation procedure
Several estimation methods were considered to estimate the imputed values of missing rainfall data. The normal ratio (NR) method has been chosen due to its formulation that can be implemented through the MI approach. Besides that, NR is selected due to its simplicity and efficiency in estimating missing rainfall values, as observed in several studies such as Mair and Fares (2010), Yozgatligil et al. (2013), and Radi et al. (2015). Several modified versions of NR methods proposed by Burhanuddin et al. (2017a) were considered in this study. They have robustified the existing ONR method by considering the trimmed mean and geometric median to improve the accuracy of the estimation results. One more modified NR was introduced in this study, i.e. NR based on the median. The estimation methods were implemented using both MI-boot and MI-blockboot approaches discussed earlier. The performance of the methods was compared to identify the best one. The different estimation methods are detailed as follows: Old normal ratio (ONR) The existing old normal ratio (ONR) estimation method was modified to improve its performance in handling the missing rainfall data and, thereby treat the outliers to reduce their effect on the estimation results. Paulhus and Kohler (1952) were the pioneers to apply the ONR method in estimating the missing rainfall data. The method is given as follows: where t and i are the sample means of the available data at target station t and i th neighbouring station respectively; t Y is the missing data at target station t; i Y is the concurrently observed data at the i th neighbouring station; N is the number of surrounding stations.

Modified normal ratio based on the trimmed mean (NRTR)
Trimmed mean is used to replace the arithmetic mean in the modified NR method. A 5% trimming percentage has been considered in this study. The modification was made to Eq. (1) and expressed as follows: where t trim and i trim are the sample trimmed means of the available data at target station t and i th neighbouring station respectively Modified normal ratio based on the median (NRMED) The mean or median is a simple way to impute the missing values. However, since the mean is very sensitive to the existence of outliers, the median is suggested to assure robustness (Fukuda and Rosta, 2005). Median imputation is preferable when the distribution of the underlying variable is not symmetric but rather skewed. Saeed et al. (2016) and Khamkong et al. (2017) applied the median to impute the missing values in rainfall dataset for their research studies. However, in this study, the median was considered as the weighting factor that replaces arithmetic mean in Eq. (1): where t med and i med are the sample medians of the available data at target station t and i th neighbouring station respectively Modified normal ratio based on the geometric median (NRGMED) The final modified NR method considered geometric mean (Gmed) in an effort to improve its performance. The Gmed is defined as the data minimizing the sum of distances to the sample dataset. The method is defined as follows (Das and Imon, 2014): Finally, Gmed is obtained by the exponential of the median of logarithm values.

Gmed exp( ) A
For this study, the Gmed is introduced in NR replacing the arithmetic mean. The formulation of the NRGMED method is as follows: where Gmed t and Gmed i are the medians of the available data at target station t and i th neighbouring station respectively The application of the estimation methods through both MI-boot and MI-blockboot approaches is implemented based on the procedure given below. All the existing and modified NR methods were considered to be implemented and tested. The ONR method was chosen to be discussed as an example.
i. The new bootstrapped samples with created artificial missing values were utilized. The procedure of introducing missing data in the dataset will be explained in the following Application section. ii. The ratio means of each sample of the target station and the dataset of its nearby stations are considered as the weighting factors for the estimation methods (an example of the ONR methodrefer to Eq. (1)). iii. Five weighting factors are obtained to produce five different estimated values for each individual missing value. iv. Each estimated value will be used to impute the missing values in the original rainfall dataset, thus, producing five complete datasets with different imputed values and the same observed values. v. The complete imputed datasets are analysed and the results of the analysis are pooled to produce a single result of the analysis.

Performance comparison
The performance of estimation methods in terms of error measurement is determined to evaluate the accuracy of the estimation results produced. Five elements of performance criteria were considered for this purpose, i.e. MSE, RMSE, MAE, and SIndex. The procedure of evaluation is executed as follows:

Results and discussion
3.1. Summary of rainfall pattern Fig. 4 shows the matrix plot of the pattern of average monsoonal rainfall amount for all the rainfall stations considered in this study. The data on the amount of average rainfall for a period of 40 years for four monsoon periods i.e. southwest, inter-monsoon 1, northeast, and inter-monsoon 2 is comprehended from the figure. Southwest monsoon occurred on four consecutive months starting from May to August while northeast monsoon prevails from November to February. Inter-monsoon 1 represented the transition of southwest monsoon while inter-monsoon 2 is the transition of the northeast monsoon.
It can be seen that the rainfall event in the Southern region was exposed to the northeast monsoon. The amount of rainfall received during these monsoons was rather significantly high, especially for Bt. 42 Jln Kluang/ Mersing and Simpang Mawai -Kuala Sedili stations. This region received a rather high amount of rainfall during this monsoon i.e. up to 1050 mm. Meanwhile, for the West region, the rainfall event was not influenced by monsoon seasons. The pattern did not exhibit an obvious difference with the amount of rainfall received by each monsoon. The rainfall received by the stations in this region was distributed evenly throughout the year.   5 displays the matrix plot of the average monthly rainfall data pattern for all of the stations in both regions. It was observed that the level of average monthly rainfall received by the Southern region appears to be high at the end of the year (December). The amount of rainfall received by this region ranges between 97 mm to 435 mm with the highest at Bt. 42, Jln Kluang/ Mersing. The month of December is usually associated with the northeast monsoon season in which the stations will receive a higher amount of rainfall compared to the other months. For the West region, the average amount of monthly rainfall received by each station is observed to be evenly distributed throughout the year. The rainfall amount received by this region ranged between 24 mm and 310 mm. Therefore, it can be concluded that the rainfall event in the Western region was less exposed to the monsoon season.

Performance evaluation of imputation approaches
The errors between the estimation results and the observed values were obtained to assess the performance of the NR methods through their implementation of both MI-boot and MI-blockboot. The most appropriate estimation methods are determined based on the least value of MSE, RMSE, and MAE and the highest value of SIndex. The comparison of the performance of the association of NR methods and MI approaches for Dataset 1 and Dataset 2 to assess their capability in both situations, i.e. without and with the presence of created outliers in the dataset, respectively was carried out. The performance evaluations have been comprehensively explained below for each estimation method and rainfall station. Fig. 6 shows the performance comparison of NR methods implemented through both MI-boot and MI-blockboot on Dataset 1 for Johor Bahru (a) and Lalang Sg Lui stations (b). The results illustrated correspond to the performance with a data missingness of 5% based on MSE, RMSE, MAE, and SIndex. The ONR, NRTR, and NRMED methods implemented through MI-blockboot are observed to provide more accurate estimation results compared to the MI-boot for both target stations. Block bootstrap has presented its benefits when dealing with time series data by providing more precise information on the data characteristics. The advantage of the combination of robustness and blocking elements has successfully improved the accuracy of the estimation results. The adoption of trimmed mean and median as a robust estimator can be clearly witnessed for the application at both stations.
NRGMED method seemed to perform even better when implemented using MI-boot compared to MI-blockboot at Lalang Sg Lui station. This method was expected to perform better when incorporated with block bootstrap similar to the other methods. Therefore, in order to identify the causes of this situation, this study was going to execute all NR methods on the dataset that was affected by the presence of outliers (i.e. Dataset 2). The purpose of employing all the methods was to prove the capability of robust methods in estimating more accurate results. The efficiency of the methods was expected to be revealed when dealing with outliers. The strength of this method in producing more accurate estimation results can be more highlighted when dealing with dataset 2.
The consistency of the performance (MSE, RMSE, MAE, and SIndex) of each method for different levels of missingness (%) was evaluated and plotted in Fig. 7 for Johor Bahru and and Fig. 8 for Lalang Sg Lui stations. The performance of all methods for both MI-boot and MI-blockboot for Johor Bahru station slightly decreased with an increase in the level of missing data percentages. The comparison between the MI approaches indicate more accurate estimation results were produced by the controlled sampling MI approach for all levels of missingness (%).
Meanwhile, it can be noted that the performance of the NR methods for Lalang Sg Lui station was rather consistent as they were not really affected by the increasing number of missing values in the dataset. The MI-blockboot was found to perform consistently better than the other MI approach, except for NRGMED method. The NRGMED incorporated MI-boot approach provided more accurate estimation results irrespective of the level of missingness (%). The robustness of an estimation method and the imputation approach is measured by its capability to be sturdy and insensitive to the presence of outliers in a dataset and maintain its performance by producing accurate results (Khamkong et al., 2017). The performance of NR methods implemented using both of the MI approaches was tested for different levels of outliers (%). Three levels of created outliers, i.e. 5%, 10%, and 15% were considered in Dataset 2 (as discussed in Section 4). Table 3 and Table 4 present the performance comparison of different NR methods implemented through MI-boot and MI-blockboot on Dataset 2 for Johor Bahru and Lalang Sg Lui stations, respectively. In order to detail examine the effects of outliers, the discussion will be firstly focused on the results of imputation at 5% missingness. Overall, it is shown that the consideration of the various number of outliers in datasets has tremendously affected the estimation results. Their performance declined with increasing percentage of outliers and the effect was significant on MI-boot's method results, especially the NRGMED. The method exhibited highly sensitive to the presence of outliers; thus, its performance was seriously influenced by the increasing number of outliers, for both Johor Bahru and Lalang Sg Lui stations. Based on the four elements of criteria, MI-blockboot's methods produced a more accurate estimation result compared to the MI-boot, in spite of increasing outliers' percentage (refer Table 3 and Table 4, the values in bold represent more accurate results produced). Table 3 demonstrates that all NR methods implemented using MI-blockboot performed better compared to the MI-boot for all levels of outliers' percentage, particularly the NRMED and NRGMED methods. Those methods were revealed to remain stable and powerful even under in the influence of outliers. The strength of the robust NRGMED associated with the controlled sampling (block bootstrap) was very obvious compared to its combination with conventional sampling (general bootstrap) as it continued to produce robust estimation results with increasing outliers in dataset. Meanwhile, the ONR and NRTR methods were seen as a sensitive and less stable option owing to their vulnerability to be affected with an increase in the outliers for both MI approaches, specifically at 15%. Therefore, the MIblockboot's NRMED and NRGMED were found to be the most appropriate methods for Johor Bahru station.
Meanwhile, based on Table 4, through the MI-blockboot implementation, all methods produced more accurate estimation results for all levels of outliers' percent compared to the MI-boot. The NRGMED was identified to perform better than other methods through MI-blockboot. The method produced robust estimation results even when dealing with the dataset that was affected by a large number of outliers. Besides that, through the implementation of both MI approaches, ONR method has performed well over the NRTR and NRMED for higher level of outliers (10% and 15%). The capability of the methods was rather affected by the presence of large number of outliers, especially the NRMED. Thus, the most appropriate methods for imputing the missing data in Lalang Sg Lui station were the NRGMED and the ONR, as they exhibited a higher tendency to be robust.
The consistency of the performance for each method for various levels of missingness (%) was also evaluated (refer to Table 3 and Table 4). The evaluation was also made for different levels of outliers. As seen in the nature of performance for Dataset 1, all the methods performed with an acceptable level of consistency with slight increment and decrement for increasing levels of missingness (%). The performance of all the estimation methods implemented using both MI-boot and MI-blockboot were considered consistent as they were not really affected by the increasing number of missing values in the dataset.  From the previous discussion, it can be concluded that the association of the robust NR methods with the controlled sampling approach (block bootstrap) has provided a more reliable and efficient approach in dealing with time series datasets, especially the datasets of poor-quality (contains large number of outliers). The application of robust estimation methods has successfully reduced the effect of outliers on the estimated values to produced accurate estimation results. In comparison to the MIblockboot approach, MI-boot has been found less efficient and more sensitive in dealing with time series dataset, especially for the poor-quality dataset. Therefore, it can be concluded that the controlled sampling approach was able to provide information with more accuracy on the original structure of rainfall time series. Table 5 presents the most appropriate NR method to be implemented through the MI approach for each target station and type of dataset. The MI-blockboot's ONR was recommended as the best method to be executed for Dataset 1 (dataset that was not affected by outliers) while the MI-blockboot's NRGMED was suggested for Dataset 2 (dataset that was contaminated with outliers) for all levels of outliers. The suggested methods are proposed for imputation of missing rainfall data for both Johor Bahru and Lalang Sg Lui stations.

Conclusions
The original structure of a dataset has to be inevitably considered as it provides crucial information on its characteristics, especially for time series data. Thus, the consideration of block bootstrap (controlled sampling) in the development of multiple imputation was a special algorithm created for rainfall time series. The block bootstrap ensured the original rainfall time series structure was preserved within each monsoon block and consequently improved the accuracy of the estimation results. These findings have successfully improved the limitation of the existing algorithm (i.e. general bootstrap) in Amelia which does not consider the structure of time series data, especially when dealing with the dataset that contains outliers (Dataset 2).
The association with the robust estimation methods has revealed its strength in this study. The advantage of MI-blockboot was evidenced over the MI-boot when dealing with the data that contains outliers. The involvement of robust NR methods, i.e. NRTR, NRMED, and NRGMED has increased its robustification to produce more accurate results. In conclusion, the introduction of controlled sampling in replacing the conventional sampling has succeeded towards improving the performance of the MI approach. Thus, the MI-blockboot was suggested for environmentalist as an alternative approach for better estimation of environmental datasets, particularly for rainfall data characterized by seasonal patterns.