Bocci, Chiara

Rocco, Emilia

English

Florence Research

An Imputation Method For Missing Covariate DataBocci, ChiaraUniversity of Florence, Department of StatisticsViale Morgagni, 5950134 Firenze, ItalyE-mail: bocci@ds.uni.itRocco, EmiliaUniversity of Florence, Department of StatisticsViale Morgagni, 5950134 Firenze, ItalyE-mail: rocco@ds.uni.itThis paper deals with the matter of applying a geoadditive model to produce estimates forsome geographical domains in the absence of point referenced geographical data. Geoadditive modelsintroduced by Kammann and Wand (2003), allow to analyze the spatial distribution of the studyvariable while accounting for possible linear or non-linear covariate eects by merging an additivemodel (Hastie and Tibshirani, 1990) and a kriging model (Cressie, 1993) and by expressing both as alinear mixed model. Therefore, when data are spatially located and explicit consideration is given tothe possible importance of their spatial distribution in the analysis or in the interpretation of results,geoadditive models represent a powerful geostatistical methodology. However, their implementationneeds the statistical units to be referenced at point locations and if we use them to produce model-basedestimates of a parameter of interest for some geographical domains, the spatial location is requiredfor all the population units. Often we don't know the exact location of all the population units,especially when socio-economic data are involved. Typically, we know the coordinates for sampledunits (which could be specically collected for the analysis), but we don't know the exact location ofall the non-sampled population units. For the non-sampled units we know just the areas to whichthey belong like census districts, blocks, municipalities, etc. In such situation, the classic approach isto locate all the units belonging to the same area by the coordinates (latitude and longitude) of thegeographical centre or centroid of the area. This is obviously an approximation, induced by nothingbut a geometrical property, and its eect on the estimates can be strong, depending on the levelof nonlinearity in the spatial pattern and on the area dimension. In this paper we propose to llthe holes in the geographical information following a stochastic imputation approach instead of theclassic deterministic one with the centroids. In particular we suggest to treat the lack of geographicalinformation imposing a distribution for the locations inside each area. This is realized through ahierarchical Bayesian formulation of the geoadditive model in which a prior distribution on the spatialcoordinates is dened. The performance of our imputation approach is evaluated through variousMarkov Chain Monte Carlo (MCMC) experiments.Methodological Framework and Basic AssumptionLet ti, 1  i  n, be a linear predictor of yi at spatial location si, s 2 R2. A geoadditive modelfor such data can be formulated as(1) yi = + tit + h(si) + "i; "i  N(0; 2");where h is an unspecied bivariate smooth functions. Representing h() with a low-rank thin plateInt. Statistical Inst.:  Proc. 58th World Statistical Congress, 2011, Dublin (Session CPS045) p.4996spline with K knotsh(s) = 0s + sTs +KXk=1ukbtps(s;k)the model (1) can be written as a mixed model (Kammann and Wand, 2003)(2) y = X + Zu+ ";withE"u"#="00#; Cov"u"#="2sIK 00 2"In#:whereX =1; ti; sTi1in ; =0; t;Ts; 0 = + 0su = [u1; :::; uK ] ;and Z is the matrix containing the spline basis functions, that isZ = [btps(si;k)]1in;1kK = [C (si   k)]1in;1kK  [C (h   k)] 1=21h;kK ;where C(v) = kvk2 log kvk and 1; :::;K are the spline knots locations.The amount of smoothing for the spline component of the model can be quantied through thevariance components ratio 2"=2s .The addition of others explicative variables is straightforward: smoothing components are addedin the random eects term Z, while linear components can be incorporated as xed eects in the Xterm. Moreover, the mixed model structure provides a unied and modular framework that allows toeasily extend the model to include various kind of generalization and evolution.Now, suppose to have a population of N units divided in Q regions, and to be interested inestimating the regional mean of a study variable y. We take a sample of n units from which we collectthe response variable y, the location s and, possibly, some other covariates (that are known withouterror for all the population units). To obtain the regional mean, we want to apply a model-basedmean estimator based on (1):^yq =1Nq24Xi2Sqyi +Xi2Rqxi^ + ziu^35 =(3)=1Nq24Xi2Sqyi +Xi2Rq ^0 + ^tti + sTiq^s +KXk=1u^kbtps(siq;k)!35 ;where Nq is the total number of units in region q, q = 1; :::; Q, and Sq and Rq indicate respectivelythe indexes of the sampled units and of the non-sampled units belonging to region q.We obtain the estimated parameters from the sampled units, but we cannot use directly (3) aswe don't know s for the not-sample units. In the classic approach the siq values are replaced withthe region centroid cq, that is a constant for all the units in region q. Here we suggest to imputethe si values through a stochastic Bayesian imputation procedure. Thus we adopt a hierarchicalBayesian formulation of model (1) (Ruppert et al, 2003) with a prior distribution fs(q) for si insideeach region q and then use the joint posterior distribution of all parameters given the data as theInt. Statistical Inst.:  Proc. 58th World Statistical Congress, 2011, Dublin (Session CPS045) p.4997basis of inference. Thus, under stochastic imputation, our complete hierarchical Bayesian formulation(following specications of Crainiceanu at al, 2003) is(4)yij;u; 2" ind N 0 + tti + Ts s+KXk=1ukbtps(s;k); 2"!;uj2sN(0; 2sIK);sijqfs(q);with non-informative priors for q, not specied as depending on the choice of fs, and for , 2s , 2"8<:0; t; sind N(0; 108) 2s ;  2"ind Gamma(10 8; 10 8):The parametrization of the Gamma(a,b) distribution implies that the parameter has mean a=b = 1and variance a=b2 = 108. Moreover, it should be noticed that we parameterize the inverse of thevariance, that is the precision parameter.We should note that the choice of the prior distribution fs(q) identies a specic imputationprocess. Anyway this choice may take into account the a priori knowledge about the spatial distri-bution of the studied phenomenon and the spatial distribution of the sampled points (that is theirspatial representativeness). If the true function fs is known, it is obviously "opportune" to use it, butthis assumption is against our work hypothesis. If no information on the distribution of the true pointlocations is available but the sample can be assumed representative of it, a non-informative fs whichlearns from the sample is a "good" choice. In particular if the spatial pattern is not too complex(like multi-modal or clustered distributions) the Beta distribution can be used. When the spatialpattern is more complex a more exible (i.e. mixture model) prior should be used. However the useof a non-informative prior works well as long as the spatial distribution of the sampled points reectsthe spatial distribution of the population points. When this condition is not satised, because thesampling design is not self-weighting and the inclusion probabilities depend on the spatial locations,a possible prior is a categorical distribution with each category corresponding to the coordinates ofa sampled point and associated probability proportional to the corresponding sampling weights. Thesame solution could be adopted if the spatial non-representativeness of the sample does not dependon the sampling design but on non sampling errors as long as the survey weights can be estimated.When the relation between the spatial distribution of the sampled points and the spatial distributionof the population is unknown, the problem can only be solved making assumptions on this relation.Finally, we note that if for each region fs(q) is a probability mass function that assumes value 1 whens = centroid and otherwise 0, then our formulation corresponds to the centroid imputation approach.MCMC ExperimentsIn order to evaluate the performance of our approach with respect to the centroids classicapproach, various MCMC experiments are implemented under dierent scenarios. All the analysis areimplemented using the OpenBUGS Bayesian inference package. We access OpenBUGS using the packageBRugs in the R computing environment.All MCMC scenarios are characterized by the following setting: The study variable is simulated by the modelyi = + xxi + f(si) + "iwhere "i  N(0; 2"), " = 0:2,  = 10, x = 0:4, x  Ber(0:5) is a dummy variable known for thewhole population, s represents the spatial location that is generated by a dierent spatial pointInt. Statistical Inst.:  Proc. 58th World Statistical Congress, 2011, Dublin (Session CPS045) p.4998process in each scenario and function f(s) is obtained as a bivariate normal mixture density andis represented in Figure 1(a). The population consisting of N = 3000 units is located in the unit squared O = [0; 1] [0; 1] which is divided in Q = 9 rectangular regions that can be represented by their vertices[(l1q;m1q); (l2q;m1q); (l2q;m2q); (l1q;m2q)]. The regions are obtained using a random binarysplitting procedure. f(s) is modeled considering a penalized thin plate spline function with Ks = 64 knots selectedon a regular grid on space. We choose this type of splines since it tends to have good numericalproperties and, as pointed out by Crainiceanu et al.(2005, p.2), the posterior correlation ofparameters for the thin-plate splines is much smaller than for other basis, which greatly improvesmixing. the MCMC analysis is implemented with a burn-in period of 15000 iterations and then weretain 5000 iterations, thinned by a factor of 5, resulting in a sample of size h = 1000 retainedfor inference. the geoadditive model (4) is tted using a sample of n = 500 units and in order to take intoaccount not only the model variability but also the variability due to the sampling design eachMCMC experiment is repeated 100 times.Each scenario diers from the others for the spatial point process used to generate s (an ho-mogenous Poisson process, an inhomogeneous Poisson process and a cluster Poisson process) and forthe sampling design (proportional stratied simple random sampling (PSRT) with strata correspond-ing to the regions and unequal probability (UP) design). Figures 1(b), 1(c) and 1(d) show the threedierent spatial patterns. Figure 2 shows the posterior density of the regional mean estimator underfour dierent scenarios: (A) an homogeneous Poisson process for s and a PSTR sampling design; (B)an inhomogeneous Poisson process for s and a PSTR sampling design; (C) a cluster Poisson processfor s and a PSTR sampling design; and (D) an homogenous Poisson process for s and an UP samplingdesign.For each scenario the regional mean estimator is evaluated under three dierent types of im-putation, that is considering three dierent a priori distributions for s: uniform distribution, betadistribution and mass function on the centroid. In addition specically for the scenario D we considerthe imputation using a categorical distribution and the mean estimator (3) is modied in order toweight the observed y values with their corresponding sampling weights. For lack of space the Figure2 shows the results only for 3 of the 9 geographic domains in which the study area is partitioned.Concluding remarksIn the last years the use of geostatistical techniques to produce model-based estimates of aparameter of interest for some geographical domains is grown. Their use however is not alwaysstraightforward as it needs for all the population units to be referenced at point location, but thisrequirement in not so easy to be accomplished. In this paper we suggest a solution to this problemthat propose a hierarchical Bayesian formulation of a geoadditive model in which a prior distributionfor the spatial coordinates is dened. The missing spatial coordinates are then extracted from theirposterior distribution, obtained by MCMC simulation.Observing Figure 2 it is straightforward to note that, when the sampling design is self-weighting,if the imputation distribution corresponds to the population spatial distribution, the stochastic impu-tation approach produces better estimates than the classic centroid approach. This is the case of theInt. Statistical Inst.:  Proc. 58th World Statistical Congress, 2011, Dublin (Session CPS045) p.49990.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0s[1]s[2]*********123456789(a)0.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0s[1]s[2]123456789(b)0.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0s[1]s[2]123456789(c)0.0 0.2 0.4 0.6 0.8 1.00.00.20.40.60.81.0s[1]s[2]123456789(d)Figure 1: Spatial distributions of the population units: (a) homogeneous Poisson process, (b) inhomogeneousPoisson process, (c) inhomogeneous Poisson process on each region, (d) independent bivariate Beta distributionon each region.Uniform approach in scenario A. Also the Beta imputation approach works well in scenario A, due tothe fact that the true spatial distribution in each region is a special case of the bivariate Beta distri-bution, but it produces less precise estimates than the Uniform imputation since the Beta parametersneed to be estimated in the tting process. When none of the imputation models corresponds to thepopulation spatial distribution, the Beta approach still presents a good performance. This dependson the fact that the Beta distribution has the advantage of modeling dierent shapes depending onthe parameters value. In our approach these parameters are estimated directly in the MCMC pro-cess exploiting the spatial distribution of the sampled units and producing a posterior bivariate Betadistribution that is as similar as possible to the sample spatial distribution. Obviously, the good per-formance of this approach relies on the representativeness of the sample, as seen in scenario D. Undersuch scenario, we obtain a good performance with the Categorical imputation approach that considerexplicitly the sampling weights.Finally when the spatial pattern is more complex the use of a more exible prior may pro-duce better results then the Beta imputation approach. Research in the use of a mixture of Betadistributions would certainly be our next step.Int. Statistical Inst.:  Proc. 58th World Statistical Congress, 2011, Dublin (Session CPS045) p.500010.0 10.2 10.4 10.6 10.8051015region 3 meanyBarposterior density10.6 10.8 11.0 11.20246810region 4 meanyBarposterior density10.6 10.8 11.0 11.202468region 5 meanyBarposterior density(a) Scenario A10.0 10.4 10.80246810region 3 meanyBarposterior density10.7 10.8 10.9 11.0 11.1 11.2051015region 4 meanyBarposterior density10.4 10.6 10.8 11.0 11.202468region 5 meanyBarposterior density(b) Scenario B9.8 10.0 10.2 10.4 10.6051015region 3 meanyBarposterior density10.6 10.8 11.0 11.202468region 4 meanyBarposterior density10.4 10.6 10.8 11.0 11.202468region 5 meanyBarposterior density(c) Scenario C10.0 10.2 10.4 10.60246812region 3 meanyBarposterior density10.6 10.8 11.0 11.20246810region 4 meanyBarposterior density10.4 10.6 10.8 11.0 11.202468region 5 meanyBarposterior density(d) Scenario DFigure 2: Posterior density of the regional model-based mean estimator under the four scenarios and for thefour imputation approaches: Centroid (green line), Uniform (red line), Beta (blue line) and Categorical (purpleline). The vertical lines indicate the true mean values.REFERENCESCrainiceanu, C., Ruppert, D., Wand, M.P. (2005). Bayesian analysis for penalized spline regression usingWinBUGS. Journal of Statistical Software 14.Cressie, N. (1993). Statistics for Spatial Data (revised edition). Waley, New York.Hastie, T.J., Tibshirani, R. (1990). Generalized Additive Models. Chapman & Hall, London.Kammann, E.E., Wand, M.P. (2003). Geoadditive Models. Applied Statistics 52, 1-18.Ruppert, D., Wand, M.P., Carroll, R.J. (2003). Semiparametric Regression. Cambridge University Press,Cambridge.Int. Statistical Inst.:  Proc. 58th World Statistical Congress, 2011, Dublin (Session CPS045) p.5001

An Imputation Method For Missing Covariate Data

https://flore.unifi.it/bitstream/2158/520274/2/Bocci_Rocco_isi2011.pdf

An Imputation Method For Missing Covariate Data

Abstract

Similar works

Full text

Available Versions

Florence Research