Estimation of missing values in aggregate level spatial data

Spatial regression models^, namely spatial lag model (SLM), spatial error model (SEM), spatial Durbin model (SDM), spatial Durbin error model (SDEM) and spatial lag of X model (SLX) were fitted in order to impute the missing values using the stochastic regression imputation method. The description of spatial regression models under comparison is given below. Let us consider the general nesting spatial model,



where $y$ is an $N\times 1$ vector of values of the dependent variable, W is an $N\times N$ dimensional neighbors weights matrix, where all elements $w_{ij}>0$ for all neighbouring units $i$ and $j$ $(i\ne j)$, and zero otherwise. $X$ is an $N\times K$ matrix of $K$ independent variables, $\beta$ is a $K\times 1$ vector of regression coefficients, $\mathrm{\rho }$ represents the autoregressive scalar parameter in the dependent variable, $\theta$ is $K\times 1$ vector of spatial spill-over parameters, $u$ represents an $N\times 1$ vector of spatially autocorrelated disturbances, $\mathrm{\lambda }$ represents the autoregressive scalar parameter in the disturbances. The parameter vector $\mathrm{\gamma }$ specifies the correlation between $X$ and the disturbance term vector $u$, ${\delta }_k$ represents the autocorrelation in the $k^{th}$ independent variable. $v_k$ and $\mathrm{\epsilon }$ are independent and randomly distributed disturbances following $\text{N }(0,{\sigma }_v^2)$ and $\text{N }(0,{\sigma }_{\epsilon }^2)$, and $x_k$ is the $k$th column-vector of $X$ for $k=1,\dots ,K$ independent variables.

The reduced form of this model can be written as,

When $\theta =0$ and $\lambda =0$ in (1), the model takes the form of spatial lag model (SLM). The SLM assumes that the dependent variable in one unit is influenced by the spatially weighted dependent variable of neighbouring units. Spatial error model (SEM) is a another form of spatial model which is formed when $\rho =0$ and $\theta =0$. In this specification it is assumed that the spatial association among the units is produced due to the unobserved features, which are either random or consist of a spatial pattern, and are independent of the explanatory variables included in the model. Spatial Durbin model (SDM) specification is formed when the autocorrelation is contributed by both dependent variable and independent variables, i.e. when $\lambda =0$. A specification which directly models the spatial spill-over effects by including the spatially lagged independent variables into the regression equation is known as the spatial lag of X (SLX) model, formed when $\rho =0$ and $\lambda =0$. The model which combines the specifications of SEM and SLX is known as the spatial Durbin error model (SDEM) which is formed when $\rho =0$.

One of the most simple and commonly used imputation technique in the non-spatial context is stochastic regression imputation. This method uses the ordinary least squares (OLS) regression to predict missing values and adds a normally distributed residual term to each predicted value. It restores loss in variability and bias associated with regression imputation (Little and Rubin 2019). The stochastic regression imputation is based on complete cases and is given as,

where $z_i$ is a random draw from $N(0,{\sigma }^2)$ and ${\sigma }^2$ is assume to be constant for given value of $x_j$. The variable with missing values $Y_k$ is the variable of interest. This variable is considered as the dependent variable of the stochastic regression imputation model and the auxiliary variables $(x_1,\dots ,x_{k-1})$ used to estimate the missing values are regarded as the independent variables.

In the present study, each missing value in dependent variable was estimated by stochastic regression imputation based on above mentioned five different spatial models. The performance of these five spatial models in accurately estimating the missing values were compared using root mean square error. For comparison purpose the stochastic regression based on OLS was also performed for estimating the missing values. Both simulated data and a real data were used to assess the performance of these six models (SLM, SEM, SDM, SDEM, SLX and OLS) in accurately estimating the missing values in aggregate level spatial data. Simulation was done to study the performance of these modes at varying percentage of missing values as well as for varied values of autocorrelation of dependent and independent variables in the model. Finally, the performance of the models were assessed using a real data from Annual Health Survey 2012–13 published by the Ministry of Home Affairs, Government of India.

Data simulation was performed using spdep package in R software to generate variables and impute missing data at aggregate level using stochastic regression imputation based on ordinary least squares (OLS), SLM, SEM, SDM, SDEM and SLX. The geographical space chosen was the map of India consisting of nine states namely Bihar, Chhattisgarh, Jharkhand, Madhya Pradesh, Odisha, Rajasthan, Uttarakhand, Uttar Pradesh and Assam aggregated at the district level. From the shape file of India with 284 districts, neighbourhood matrix W was constructed based on queen contiguity weights matrix, which defines the neighbors as those areas that share boundary points as well as some portion of their boundary. Simulated data of size 284 districts included three independent variables for which the regression coefficient vector was fixed arbitrarily at $\mathrm{\beta }={\left(\begin{array}{ccc}-3& 3.5& -3\end{array}\right)}^T$. The disturbance parameters were generated as an independent and randomly distributed normal variable with mean zero and fixed at ${\sigma }_v^2,{\sigma }_{\epsilon }^2=1$. The spatial spill-over parameter vector was fixed at $\mathrm{\theta }={\left(\begin{array}{ccc}-1& 1& -1\end{array}\right)}^T$. The autocorrelation in disturbance vector $u$ and the omitted variable bias were fixed at $\lambda =\left\{0.05\right\}$ and $\mathrm{\gamma }={\left(\begin{array}{ccc}0& 0& 0\end{array}\right)}^T$ respectively. To generate dependent variable with low and high autocorrelation, the scalar parameter $\mathrm{\rho }$ was fixed at two values namely $\mathrm{\rho }=0.2$ and $\mathrm{\rho }=0.8$ respectively. Similarly, the autocorrelation in independent variables was assumed to be $\delta =0.1$ and $\delta =0.6$ to generate independent variables with low and high autocorrelation respectively. The variable of interest (dependent variable) was finally simulated following the general nesting spatial model given in equation (3). The simulation resulted in four datasets each of size 284 representing four different situations as given below.

• a)
  High autocorrelation in both dependent and independent variables
• b)
  Low autocorrelation in both dependent and independent variables
• c)
  High autocorrelation in dependent variable and low autocorrelation in independent variables
• d)
  Low autocorrelation in dependent variable and high autocorrelation in independent variables

The autocorrelation and Pearson's correlation coefficient values of the simulated dependent and independent variables are given in Table 1.

The data was simulated 1000 times for each of the above four scenarios of levels of autocorrelation among the variables. In each simulated dataset a fixed proportion of values in dependent variable was assumed to be missing according to missing at random (MAR) mechanism. In MAR mechanism, probability of missing values in the dependent variable (Y) was based on the values of other variable(s). To generate missing as per MAR mechanism in dependent variable, a fixed proportion of locations were randomly selected among those locations where the values of the independent variable $X_1$ were very low (less than the first quartile value of$X_1$). In these selected locations the values of dependent variable were assumed to be missing. This was repeated for various proportions of missing, 5%, 10%, 15%, 20% and 25%in the dependent variable in each simulated dataset. For each proportion of missing and for each level of autocorrelation in variables, the missing values were estimated using stochastic regression imputation with ordinary least squares method (OLS), SLM, SEM, SDM, SDEM and SLX.

The performance of the above models in accurately estimating the missing values under each scenario were compared based on the computed root mean squared error (RMSE) values. For $i=1tow$ where$i$ are the regions consisting of missing values and $w$ are the total number of regions with values to be imputed, the RMSE for the variable Y is given as,

where, ${\hat{\mathit{Y}}}_{\mathit{ij}}$ is the imputed value for the variable Y in region $i$, for a given simulated dataset (j = 1, 2, …,1000), for any given model and any fixed proportion of missing in the dataset. ${\mathit{Y}}_{\mathit{ij}}$is the observed value of the variable Y in region $i$ and jth simulated dataset.

The mean of the RMSE (mRMSE) and standard deviation of RMSE (sdRMSE) for a given model and fixed proportion of missing were calculated as,

The mRMSE was computed for the imputations performed using each of the five spatial regression models and OLS model, on the data with 5%–25% proportions of missing values in the dependent variable. The mRMSE values were compared between the models and the model with lowest mRMSE in each scenario was considered the best in estimating the missing values in the dependent variable.

In reality, missing values can occur for the regions where data on certain events/diseases cannot be collected due to constraints such as safety, logistics etc. To reflect such situations, the proposed imputation techniques were applied on data obtained from AHS for 284 districts of India. The variables have been chosen in a manner to reflect the relationships in two scenarios, (i) high autocorrelation in dependent and independent variables (high-high) and (ii) low and high autocorrelation in dependent and independent variables (low-high) respectively. For the first scenario, the dependent variable considered was total fertility rate of a district (TFR) and the independent variables considered were percentage of mothers who received any antenatal check-up and unmet need for spacing in a district. In the case of second scenario namely low-low autocorrelation, the dependent variable considered was percentage of women aged 15–19 years who were already mothers/pregnant at the time of survey while independent variables total fertility rate and percentage of female literacy rate in a district. In real data analysis, overall rates of missingness of 5%, 10%, 15%, 20% and 25% of observations were imposed on the dependent variable according to MAR mechanism. Missing values were imposed on the dependent variable using the same procedure as it was done on the simulated data. The variables considered for the analysis and the values of their spatial autocorrelation and Pearson's correlation coefficient are listed in Table 2.

The stochastic regression imputation with all the six models were applied on the simulated data under two scenarios (i) for four types of levels of autocorrelation among the dependent and independent variables as mentioned in Table 1, and (ii) for various proportion of missing values in the dependent variable. Fig. 1 gives the comparison of mean RMSE values for the six models for the above mentioned two scenarios. When the autocorrelation in both dependent and independent variables were high (Fig. 1(a)), it can be noted that the imputation performed using SEM had higher mRMSE compared to that obtained using the remaining five models. The mRMSE value using SEM was even higher for the imputation performed using OLS model, which model does not account for autocorrelation. The least and similar mRMSE in the imputation were produced by SLM, SDM and SDEM. When the autocorrelation in both dependent and independent variables were low (Fig. 1(b)), the performance of imputation using all the models was more or less similar in terms of mRMSE. However, the imputations performed using SDM, SDEM and SLX found better with least mRMSE. In case of high autocorrelation in dependent variable and low in independent variables (Fig. 1(c)), the pattern of mRMSE is similar to that of situation (a), but it can be noted that the mRMSE of imputation using SEM has reduced. Imputation performed using SLM, SDM and SDEM had least mRMSE when the autocorrelation in dependent variable was low and high in independent variables (Fig. 1(d)). In this case, OLS and SEM performed in a similar manner. Overall, the simulation analysis revealed that when autocorrelation is high in the variable with missing values, stochastic regression imputation performed using spatial lag model, spatial Durbin model and spatial Durbin error model gives accurate and almost similar estimates of missing values irrespective of whether the independent variables in the regression model were spatially autocorrelated or not, and this result was consistent for the varied amount of missing values. When the autocorrelation is low in the variable with missing values, in addition to the above three models, the stochastic regression imputation using spatial lag X model was also found to be equally accurate in estimating the missing values. It can also be noted that stochastic regression imputation performed using OLS model and SEM had the highest mRMSE in all the scenarios considered in the simulation study.

The mean and standard deviations of RMSE values obtained for 1000 simulated datasets and for real data analysis in each scenario is given in supplementary material (Table 3 and Table 4).

Fig. 2 presents a line graph depicting the mRMSE obtained from imputations performed by stochastic regression using different models. Similar to the results obtained from the simulated data, imputation performed using SLM, SDM and SDEM were found to have least RMSE when the autocorrelation was high in both dependent and independent variables. This pattern remained consistent for varying proportions of missing values. The same three models performed well when the data possessed low autocorrelation in dependent variable and high autocorrelation in independent variables. It can also be noted that stochastic regression imputation performed using OLS model had highest mRMSE in all the scenarios considered.