Nonparametric imputation by data depth

1.8 MB·

by  Pavlo Mozharovskyi

Nonparametric imputation by data depth
Pavlo Mozharovskyi∗
CREST, Ensai, Universite´ Bretagne Loire
and
Julie Josse
7
´
1 Ecole Polytechnique, CMAP
0
and
2

n Franc¸ois Husson
a
J IRMAR, Applied Mathematics Unit, Agrocampus Ouest, Rennes

2
1
January 12, 2017

]
E
M
.
t Abstract
a
t
s
The presented methodology for single imputation of missing values borrows the idea from
[
data depth — a measure of centrality defined for an arbitrary point of the space with respect to
1
a probability distribution or a data cloud. This consists in iterative maximization of the depth
v
3 of each observation with missing values, and can be employed with any properly defined sta-
1 tistical depth function. On each single iteration, imputation is narrowed down to optimization
5
of quadratic, linear, or quasiconcave function being solved analytically, by linear programming,
3
0 or the Nelder-Mead method, respectively. Being able to grasp the underlying data topology, the
.
1 procedure is distribution free, allows to impute close to the data, preserves prediction possibili-
0 tiesdifferenttolocalimputationmethods(k-nearestneighbors,randomforest),andhasattractive
7
robustness and asymptotic properties under elliptical symmetry. It is shown that its particular
1
: case — when using Mahalanobis depth — has direct connection to well known treatments for
v
i multivariate normal model, such as iterated regression or regularized PCA. The methodology is
X
extended to the multiple imputation for data stemming from an elliptically symmetric distribu-
r
a tion. Simulation and real data studies positively contrast the procedure with existing popular
alternatives. ThemethodhasbeenimplementedasanR-package.
Keywords: Ellipticalsymmetry,Outliers,Tukeydepth,Zonoiddepth,Nonparametricimpu-
tation,Convexoptimization.
∗ThemajorpartofthisprojecthasbeenconductedduringthepostdocofPavloMozharovskyiatAgrocampusOuest
(Rennes)grantedbyCentreHenriLebesgueduetoprogramPIA-ANR-11-LABX-0020-01.
1
1 Introduction
Following the seminal idea of Tukey (1974), the concept of data depth has developed to a pow-
erful statistical methodology allowing description of data w.r.t. location, scale, and shape based
on a multivariate ordering. Today, it finds numerous applications in multivariate data analysis
(Liuetal.,1999),statisticalqualitycontrol(LiuandSingh,1993),classification(Jo¨rnsten,2004,
Langeetal.,2014),multivariateriskmeasurement(CascosandMolchanov,2007),robustlinear
programming (Bazovkin and Mosler, 2015), etc. Being able to exploit topological properties of
the data in a nonparametric way, statistical depth function further proves to be suited for imput-
ing missing values while being connected to the state of the art methods. We start with a brief
backgroundonimputation,followedbytheproposal.
1.1 Background on missing values
Theproblemofmissingvaluesexistssincetheearliestattemptsofexploitingdataasasourceof
knowledgeasitliesintrinsicallyintheprocessofobtaining,recording,andpreparationofthedata
itself. To exploit all the information present in the data set, a statistical method may be adapted
tomissingvalues,butthisrequiresdevelopingsuchaoneforeachestimatorandinferenceofin-
terest. Amoreuniversalwayistoimputemissingdatafirst,andthenapplythestatisticalmethod
ofinteresttothecompleteddataset(LittleandRubin,2002). Lastly,themultipleimputationhas
gainedalotofattention: foradatasetcontainingmissingvaluesanumberofcompleteddatasets
is generated reflecting uncertainty of the imputation process, which enables not only estimating
theparameterofinterestbutalsodrawinganinferenceonit(Rubin,1996).
Developmentofmanystatisticalmethodsstartedwiththenaturalnormalityassumption,and
imputation is not an exception here. For a multivariate normal distribution, single imputation
of the observations containing missing values can be performed by imputing the missing values
with conditional mean, where conditioning is w.r.t. observed values. This procedure makes
use of the expectation-maximization algorithm (Dempster et al., 1977) to estimate mean and
covariance matrix; see also Little and Rubin (2002). By that, imputation inherits sensitivity to
outliers and to (near) low-rank covariance matrix. To deal with these problems, uncountable
extensions of the EM framework have been developed. Another way is to directly assume low-
rank of the underlying covariance, which is followed by the methods based on the principal
component analysis (PCA), closely connected to the matrix completion literature, see Josse and
Husson (2012), Hastie et al. (2015). These methods aim at denoising the data and suppose its
special structure, knowledge (or estimation) of the rank of data, of the shape of outliers, and
noise to be normal. Both groups of methods impute on low-dimensional affine subspaces and
—bythat—aresensitivetoevenslightdeviationsfromnormality(orellipticityingeneral)and
ignore the geometry of data. Extension of single imputation methods to more general densities
consistsinsuchnonparametrictechniquesask-nearestneighbors(kNN)(seeTroyanskayaetal.,
2001, and references there in) and random forest (Stekhoven and Bu¨hlmann, 2012) imputation.
When properly tuned, these methods capture the locality of the data, but fail to extrapolate and
thuscanbeinappropriateunderthemissingatrandommechanism(MAR,Seamanetal.,2013),
i.e. exhibitinsomesensesuperfluouslylocalbehavior.
2
vol = 0.42061
6 6
2] 2]
5, ][, 5, ][,
0. 0. 4
< 4 <
s)) 1))
s s 2
X2.mi 2 X.mis
a( a( 0
n n
s. s.
ms(i 0 ms(i −2
u u
S S
w w
2[ro −2 1[ro −4
X X
−1 0 1 2 3 −1 0 1 2 3
X2[rowSums(is.na(X2.miss)) < 0.5, ][,1] X1[rowSums(is.na(X.miss1)) < 0.5, ][,1]
Figure1: Bivariatenormaldistributionwith30%MCAR(left)andwithMARinsecondcoordi-
nate for values > 3.5 (right); imputation using maximum zonoid depth (red), conditional mean
imputationusingEMestimates(blue),andrandomforestimputation(green).
1.2 Proposal
Currentarticlefillstheexistinggaponthelevelofsingleimputationforgeneralellipticallysym-
metric distribution families. We suggest a nonparametric framework that is able to grasp data
topology and, to some degree, even reflect finite sample deviations from ellipticity. For the pur-
pose of single imputation of missing values, we resort to the idea of statistical depth function
D(x|X)—ameasureofcentralityofanarbitrarypointx ∈ Rd w.r.t. ad-variaterandomvector
X. Given such a measure, we propose to maximize it for a point having missing values condi-
tionedonitsobservedvalues. Foreachpointcontainingmissingvalues,theprocedureisrepeated
iteratively to achieve stability of the solution. In this framework, the properties of the imputa-
tion are to great extent defined by the chosen notion of data depth function (see Section 2). A
nonparametric data depth allows for imputation close to the data geometry, still accounting for
theirglobalfeaturesduetothecenteranchoring. Thismakesadistinctionw.r.t. usingfullylocal
imputation methods like kNN or random forest imputation. In addition, data depth allows for
robustimputationbothinsenseofoutliers(i.e. disregardingoutlierswhenimputingpointscloser
to center) and distribution (i.e. not masking outliers). As not using it, depth-based imputation
avoidsproblemsconnectionwithestimationofthecovariancematrix.
We employ three depth functions: Mahalanobis depth, which is a natural extension of the
Mahalanobis distance and is here because of its connections to standard imputation frameworks
(see Section 4); zonoid depth, which imputes by the average of the maximal number of equally
weightedobservations;Tukeydepth,whichmaximizestheinfimumoftheportionofpointscon-
tained in the closed halfspace including the imputed point. For a well defined imputation, these
depth notions require two, one, and zero first moments of the underlying probability measure,
respectively. We highlight most important properties of the depth-based imputation right below
bymeansofexamples.
RegardFigure1(left),where150pointsstemmingfromabivariatenormaldistributionhav-
ingmeanµ = (1,1)(cid:48) andcovarianceΣ = (cid:0)(1,1)(cid:48),(1,4)(cid:48)(cid:1)areplotted,with30%oftheentries
1 1
havingbeenremovedinbothvariablesduetothemissingcompletelyatrandom(MCAR)mecha-
3
2] 2]
5, ][, 6 5, ][, 40
0. 0.
< 4 <
s)) s)) 20
s s
mi 2 mi
3. 4.
X X 0
na( 0 na(
s. s.
s(i s(i 20
m 2 m −
u − u
S S
w w
3[ro −4 4[ro −40
X X
−1 0 1 2 3 4 −20 −10 0 10 20
X3[rowSums(is.na(X3.miss)) < 0.5, ][,1] X4[rowSums(is.na(X4.miss)) < 0.5, ][,1]
Figure2: Mixtureofbivariatenormal(425points,15%MCAR)andCauchy(75point)samples
(left) and 1000 bivariate Cauchy distributed points with 15% MCAR (right). Imputation using
maximumTukeydepth(red)andconditionalmeanimputationusingEMestimates(blue).
nism;pointswithonemissingentryaredenotedbydottedlines. EM-basedimputationisdepicted
in blue. The imputed points are conditional means and lie exactly on the regression lines esti-
mate; their distribution versions are plotted as the two lines. Zonoid depth-based imputation,
picturedinred,reflectstheideathatthesamplemaynotnecessarilybenormal,withthisunsure-
nessmorestressedonthefringeofthedatacloud,whereimputedpointsdeviatefromconditional
meantowardstheunconditionalone. Indeed,awayfromthedatacenter,imputationbydatadepth
is closer to nonparametric imputation, which accounts for the local nature of data. On the other
hand,differenttothesemethods,depthimputationpreservespredictionability. InFigure1(right)
forthesamesamplefirstcoordinatehasbeenremovedfor(14)pointshavingvalue> 3.5inthe
second coordinate (MAR mechanism). The depth-based imputation manages to extrapolate to
predictmissingvalues,whilerandomforestimputationperformsasexpectedratherpoorly.
Beingrobustbothinsenseofoutliersandheavytailsconcernsthetwofollowingcases: First,
if data is polluted by outliers but coordinates are missing for a representative point, the outliers
should not alter the imputation. In Figure 2, left we plot — zoomed in — 500 points. 425 of
them stem from the same normal distribution as above, with 15% of entries MCAR; another
75 are outliers drawn from the Cauchy distribution with the same center and shape matrix and
having no missing values. As expected, conditional mean based on EM estimates (depicted in
blue) imputes rather randomly. Depth-based imputation using Tukey depth imputes in a robust
way,closeto(distribution)regressionlines,reflectinggeometryofdata. Second,ifmissingvalues
belong to the polluting distribution (or if the entire data generating process is heavy-tailed), the
imputationshouldreflectthisheavy-tailnessnottomaskapossibleoutlier. For1000pointsfrom
Cauchydistributionhaving15%ofvaluesMCAR,imputationbyTukeydepthandtheEM-based
one (for comparison) are shown in Figure 2, right. Depth-based respects the general ellipticity
andimputesclosetodistributionregressionlines.
Therestofthepaperisorganizedasfollows. InSection2westatesomeimportantdefinitions
concerningdatadepthandestablishthenotation. Section3enlightenstheproposedmethodology
ofthedepth-basedimputation,regardsitstheoreticalproperties,andsuggestsoptimizationtech-
niques. Section 4 is devoted to the theoretical investigation of the special case when employing
4
Mahalanobis depth bridging the proposed imputation technique with regression and PCA impu-
tation. Section 5 provides a comparative simulation and real data study. Section 6 extends the
proposedapproachtomultipleimputation. Section7gatherssomeusefulremarks.
2 Notation and background on data depth
ForadatasetinRd (n×dmatrix)X = {x ,...,x } = (X ,X ),wedenotebyX and
1 n obs miss obs
X itsobservedandmissingpart,respectively,i.e. thereexistsanindicatorn×dmatrixM
miss
suchthatX = X ·(1 −M)+X ·M inthesenseofelementwisemultiplication. For
obs n×d miss
a point x ∈ Rd, we denote miss(x) and obs(x) the sets of its coordinates containing missing,
respectively observed values, denoting |miss(x)| and |obs(x)| the corresponding cardinalities,
restrictingto|miss(x)|+|obs(x)| = d. Followingthesamelogic,wewritemiss(i)andobs(i)
forapointx ∈ X orevenjustmissandobsifnoconfusionarises.
i
Under missing completely at random (MCAR) we understand the mechanism where each
valueofX hasthesamemissingprobability.
BeinginthecoreofthetheoreticalresultsderivedinSection3ellipticaldistributionisdefined
asfollows(seeFangetal.(1990),andLiuandSingh(1993)inthedatadepthcontext).
Definition1. RandomvectorX ∈ Rd isdistributedellipticallysymmetricwithstrictlydecreas-
ingdensityorequivalentlyasE(µ ,Σ ,F)ifitisdistributedasX =D µ +RΛU withR ∈ R
X X X +
beinganonnegativerandomvariablestemmingfromF possessingstrictlydecreasingdensity,U
uniformlydistributedonSd−1,µ ∈ Rd,andΣ = ΛΛ(cid:48).
X
Below we briefly state definitions covering necessary material on data depth. Following the
pioneeringideaofTukey(1974),statisticaldatadepthfunctionisamapping
Rd×M → [0,1] : (x,P) (cid:55)→ D(x|P),
with M being a subset of M , the set of all probability measures on (Rd,B). It measures
0
closenessofxtothecenterofP. Further,foragivenpointxandarandomvectorX comingfrom
the probability distribution P ∈ M we denote the depth D(x|X), and D(x|X) its empirical
version. Zuo and Serfling (2000) developed axiomatic for the depth, which requires a proper
depth function to be affine invariant, maximal at the P’s center of symmetry, non-increasing
on a ray starting from any deepest point, and vanish in infinity. Additionally one may require
quasiconcavity, a restriction which will prove to be useful computationally below. Upper-level
set of the depth function to level α is called a depth (-trimmed) region D (X) = {x ∈ Rd :
α
D(x|X) ≥ α}. Forα ∈ [0,1],depthregionsformafamilyofnestedset-valuedstatistics,which
describesX w.r.t. location,scatter,andshape.
A number of depth notions emerged during the last decades; below we give definitions of
thosethreebeingusedformeanofimputationinthepresentarticle.
Definition2. Mahalanobis(1936)depthofx ∈ Rd w.r.t. X isdefinedas
DM(x|X) = (cid:0)1+(x−µ )(cid:48)Σ−1(x−µ )(cid:1)−1,
X X X
whereµ andΣ areanylocationandshapeestimatesofX.
X X
5
In the empirical version, µ and Σ are replaced by appropriate estimates of location and
X X
shape, which throughout the article are taken as the moment estimates µ = 1 (cid:80)n x and
X n i=1 i
Σ = 1 (cid:80)n (x −µ )(x −µ )(cid:48),respectively.
X n−1 i=1 i X i X
KoshevoyandMosler(1997)defineazonoidtrimmedregion—forα ∈ (0,1]—as
(cid:110)(cid:90) (cid:20) 1(cid:21) (cid:90) (cid:111)
Dz(X) = xg(x)dP(x) : g : Rd (cid:55)→ 0, measurableand g(x)dP(x) = 1
α α
Rd Rd
andforα = 0as
Dz(X) = conv(supp(X)) .
0
wheresupp(X)denotesthesupportofX andconv(A)denotesthesmallestconvexsetcontaining
A. Itsempiricalversioncanbedefinedas
n n
(cid:110)(cid:88) (cid:88) 1 (cid:111)
Dz(n)(X) = λ x : λ = 1, λ ≥ 0, αλ ≤ ∀ i ∈ {1,...,n} .
α i i i i i n
i=1 i=1
Definition3. Zonoiddepthofxw.r.t. X isdefinedas
(cid:26) sup{α : x ∈ Dz(X)} ifx ∈ conv(cid:0)supp(X)(cid:1),
Dz(x|X) = α
0 otherwise.
ForacomprehensivereferenceonthezonoiddepththereaderisreferredtoMosler(2002).
Zonoid depth tends to represent x as the average of maximum number of equally weighted
points,i.e. asaweightedmean,whichopensavarietyofconnectedmethodsincludingtheentire
classoftheweightedmeandepths,seeDyckerhoffandMosler(2011).
Definition4. Tukey(1974)depthofxw.r.t. X isdefinedas
DT(x|X) = inf{P(H) : H aclosedhalfspace, x ∈ H}.
In empirical version, probability is substituted by the portion of X. Exploiting solely the
datageometryandoptimizingovertheindicatorloss,Tukeydepthisfullynonparametric,highly
robust, and does not require moment assumptions on X. For more information on Tukey depth
andthecorrespondingtrimmedregionsseeDonohoandGasko(1992)andHallinetal.(2010).
3 Imputation by depth maximization
3.1 Main idea
Let’sstartbyregardingoneofthesimplestmethodstoimputemissingvaluesmakinguseofthe
following iterative regression imputation scheme: (1) initialize missing values arbitrary, using
mean imputation for instance; (2) impute missing values in one variable by the values predicted
bytheregressionmodelofthisvariablewiththerestingvariablestakenasexplanatoryones,(3)
iteratethroughvariablescontainingmissingvaluestillconvergence.
Mostcommonimputationmethodsarebasedonsimilarapproaches: Templetal.(2011)use
the same (robustified) iterative method, and Josse and Husson (2012) and Hastie et al. (2015)
use iterative (thresholded) singular value decomposition to impute the missing entries, whereas
Stekhoven and Bu¨hlmann (2012) iteratively replace missing cells with values fitted by random
6
forest. This scheme can even be traced back to the earliest solution to deal with missing values
suggested by Healy and Westmacott (1956) and is at the root of EM based algorithms (imple-
mented e.g. by means of the SWEEP operator for normal data) to estimate parameters despite
missingvalues(LittleandRubin,2002).
In the described above procedure, on each step, each point x having missing values
i
on a coordinate j is imputed with y the univariate conditional mean E[X|X =
i,j {1,...,d}\{j}
x ,µ = µ ,Σ = Σ ]. Afterconvergence,eachpointx withmissingvalueson
i,{1,...,d}\{j} X X X X i
miss(i)isimputedwiththemultivariateconditionalmean
E[X|X = x ,µ = µ ,Σ = Σ ] (1)
obs(i) i,obs(i) X X X X
=µ +Σ Σ−1 (cid:0)x −µ (cid:1).
Xmiss(i) Xmiss(i),obs(i) Xobs(i),obs(i) i,obs(i) Xobs(i)
Thelastexpressionistheclosed-formsolutionto
min dMah(z,µ |Σ )
X X
z ∈R|miss(i)|,z =x
miss(i) obs(i) obs(i)
withdMah(z,µ |Σ ) = (z−µ )(cid:48)Σ−1(z−µ )beingthesquaredMahalanobisdistancefrom
X X X X X
z to µ . Minimizing Mahalanobis distance can be seen as maximizing a centrality measure —
X
theMahalanobisdepth
max DMah(z|X).
z ∈R|miss(i)|,z =x
miss(i) obs(i) obs(i)
Wegeneralizethisprincipletothenotionofstatisticaldepthfunction.
We propose a unified framework based on the statistical data depth function. Having a sam-
pleX,imputeapointxcontainingmissingcoordinateswiththepointy maximizingdatadepth
conditioned on observed values x . This direct extension of the idea of conditional mean im-
obs
putationtodatadepthcanbeexpressedas
y = argmax D(z|X). (2)
zmiss∈R|miss|,zobs=xobs
Theuseofequation(2)islimitedtostrictlyquasiconcavecontinuousandnowherevanishing
depthnotions,whichisnotthecaseforintrinsicallynonparametricdepthsbestreflectingthedata
geometry. Asolutionto(2)cantriviallybenonunique,asfittingtothefinite-sampledatatopology
the depth can be non-continuous. In addition, the value of the depth function may become zero
immediatelybeyondtheconvexhullofthesupportofthedistribution,whichisjustconv(X)for
a finite sample. To circumvent these problems, we suggest to impute x having missing values
withy:
(cid:0) (cid:1)
y = ave argmin {(cid:107)u−v(cid:107)|v ∈ D (X)} , (3)
α∗
u∈Rd,u =x
obs obs
where
α∗ = inf (cid:8)α|D (X)∩{z|z ∈ Rd, z = x } = ∅(cid:9). (4)
α obs obs
α∈(0;1)
GivenasamplecontainingmissingdataX = (X ,X ). Startwithanarbitraryinitial-
obs miss
izationofallmissingvalues,byimputingwithunconditionalmeansay. Then,foreachobserva-
tion containing missing coordinates, impute them according to (3), and iterate. This imputation
byiterativemaximizationofdepthcanbesummarizedasAlgorithm1.
7
Algorithm1Singleimputation
1: function IMPUTE.DEPTH.SINGLE(X)
2: Y ← X
3: µ ← µˆ(obs)(X) (cid:46)Calculatemeanignoringmissingvalues
4: fori = 1 : ndo
5: ifmiss(i) (cid:54)= ∅then
6: y ← µ (cid:46)Imputewithunconditionalmean
i,miss(i) miss(i)
7: I ← 0
8: repeat (cid:46)Iterateuntilconvergenceormaximaliteration
9: I ← I +1
10: Z ← Y
11: fori = 1 : ndo
12: ifmiss(i) (cid:54)= ∅then (cid:46)Imputewithmaximumdepth
13: α∗ ← inf (cid:8)α|D (X)∩{z|z ∈ Rd, z = y } = ∅(cid:9)
α∈(0;1) α obs i,obs(i)
(cid:0) (cid:1)
14: yi ← ave argminu∈Rd,u =y {(cid:107)u−v(cid:107)|v ∈ Dα∗(X)}
obs i,obs(i)
15: untilmax |y −z | < (cid:15) or I = I
i∈{1,...,n},j∈{1,...,d} i,j i,j max
16: returnY
After the stopping criterion has been reached, one can expect that for each point initially
containing(cid:9)(cid:12)missing values, it holds xi = argmaxzobs=xobsminu∈Sd−1(cid:12)(cid:12)(cid:8)k : x(cid:48)ku ≥ z(cid:48)u, k ∈
{1,...,n} (cid:12) when employing the Tukey depth. So the imputation is performed due to the max-
imin principle based on criteria involving indicator functions, which implies robustness of the
solution. When using zonoid depth, each such x is imputed by the average of the maximum
i
number of possibly most equally weighted points. W.r.t. X, this is a weighted mean impu-
tation, which has connection to the methods of local nature such as the kNN imputation, and
allowsforgainingadditionalinsightsintodatageometrybyinspectingtheoptimalweights,con-
stitutingtheLagrangemultipliers;seeSection3.3forthedetaileddiscussion. WithMahalanobis
depth, each x with missingness is imputed by the conditional mean (1) and thus lies in the
i
single-output regression hyperplane X on X for all j ∈ miss(i) and in general
·,j ·,{1,...,d}\{j}
(cid:0) (cid:1)
inthe d−|miss(i)| -dimensionalmultiple-outputregressionsubspaceX onX ;
·,miss(i) ·,obs(i)
suchasubspaceisobtainedastheintersectionofthesingle-outputregressionhyperplanescorre-
spondingtomissingcoordinates. Amongothers,thisyieldsfurtherpropertiesandconnectionsto
covariancedeterminantandtheregularizedPCAimputationbyJosseandHusson(2012),which
isregardedindetailinSection4.
3.2 Properties
Though being of different nature, suggested depth-based framework maintains some of the de-
sirable properties of the EM-based imputation, since the imputation converges to center of the
conditionaldistributioninthesettingsspecifiedinTheorems1and2. DifferenttoEM,weavoid
thesecondandfirstmomentassumptionsduetotheusezonoidandTukeydepths. Togetherwith
Mahalanobis depth, the choice of depths becomes canonical in the sense that finiteness of first
two(Mahalanobisdepth),one(zonoiddepth),andno(Tukeydepth)momentsisrequired.
8
Theorem 1 shows that for any elliptical distribution, imputation of one point only converges
tothecenteroftheconditionaldistributionwhenconditioningontheobservedvalues.
Theorem 1 (One row consistency). Let X(n) = (X(n),x) be a sequence of data sets in Rd
obs
(n)
with d ≥ 2 consisting of a point x = (x ,x ) and X being n points sampled from
obs miss obs
an E(µ ,Σ ),f). Then for n → ∞ for Tukey, zonoid (E possesses finite 1st moment), and
X X
Mahalanobis(E possessesfinite2ndmoment)depths
y = µ +Σ Σ−1 (x −µ )
miss Xmiss Xmiss,obs Xobs,obs obs Xobs
isastationarypointofAlgorithm1.
Theorem1isillustratedinFigure3forabivariatesamplestemmingfromtheCauchydistri-
bution, whereapointisimputedusingTukeydepth. Thekerneldensityestimateoftheimputed
valuesover10000repetitionsresemblestheGaussiancurveandapproachesthetargetvaluebe-
ingthecenterofthepopulationconditionaldistribution.
Student t1, n = 100, d = 2, k = 10000 Student t1, n = 200, d = 2, k = 10000
mputed values 0.60.8 mputed values 0.60.8
mate of i 0.4 mate of i 0.4
esti esti
ensity 0.2 ensity 0.2
d d
el el
n n
Ker 0.0 Ker 0.0
0 S1tudent t21, n = 5030, d = 2, 4k = 100050 6 0 S1tudent t12, n = 10030, d = 2,4 k = 100050 6
mputed values 0.60.8 Maximum Tukey depth imputation of X[i,2] for X[i,1] = 3 mputed values 0.60.8 Maximum Tukey depth imputation of X[i,2] for X[i,1] = 3
mate of i 0.4 mate of i 0.4
esti esti
ensity 0.2 ensity 0.2
d d
el el
n n
Ker 0.0 Ker 0.0
0 1 2 3 4 5 6 0 1 2 3 4 5 6
Maximum Tukey depth imputation of X[i,2] for X[i,1] = 3 Maximum Tukey depth imputation of X[i,2] for X[i,1] = 3
Figure 3: Kernel density estimate (solid) and the best approximating Gaussian curve (dashed)
over10000repetitionsoftheimputationofasinglepointhavingonemissingcoordinate. Sample
ofsize100(top,left),200(top,right),500(bottom,left),1000(bottom,right)isdrawnfromthe
Cauchy distribution with location and scatter parameters µ and Σ from the introduction. The
1 1
populationconditionalcentergiventheobservedvalueequals3.
Theorem 2 states that if missing values constitute a portion of the sample but concentrate in
a single variable, the imputed values converge to the center of the conditional distribution when
conditioningontheobservedvalues.
9
Theorem 2 (One column consistency). Let X(n) = (X(n),X(n) ) be a sequence of data sets
obs miss
(n)
with x only in coordinate j following MCAR mechanism, sampled from an E(µ ,Σ ,F).
miss X X
Then for n → ∞ for Tukey, zonoid (E possesses finite 1st moment), and Mahalanobis (E pos-
sessesfinite2ndmoment)depthsforallpointshavingmissingvalues
y = µ +Σ Σ−1 (x −µ )
j Xj Xj,obs Xobs,obs obs Xobs
withobs = {1,...,d}\{j}isastationarypointofAlgorithm1.
The statements of both theorems are valid for n → ∞. Nevertheless, for any (cid:15) > 0 there
exists N such that for all n > N (cid:15)-neighborhood of the corresponding limit from Theorems 1
and2isastationarypointofAlgorithm1.
3.3 Optimization in a single iteration
Due to Algorithm 1, the only imputation step repeated iteratively consists in imputing point x
having missing coordinates miss(x) with y by maximizing its depth D(z|X) w.r.t. the com-
pleteddatasetX conditionedonx . Notethatthefunctionf(z ) := D(z|X,z =
obs miss(x) obs(x)
x )isquadraticforMahalanobisdepth,continuousinsideconv(X)forzonoiddepth,andstep-
obs
wise discrete there for Tukey depth, with the single solution to a convex maximization problem
inallthethreecases. ForatrivariateGaussiansample,f(z )isdepictedinFigure4
miss(x)
For Mahalanobis depth, equation (2) has closed form: y = µ +
miss Xmiss
Σ Σ−1 (x −µ ). In case Σ is singular, one can work in the linear sub-
Xmiss,obs Xobs,obs obs Xobs X
spaceofeigenvectorswithpositiveeigenvalues. DefinedthiswayMahalanobisdepthissensitive
to outliers, which can be compensated for by estimating µ and Σ in a robust way, using
X X
minimumcovariancedeterminant(MCD,seeRousseeuwandVanDriessen,1999)say.
Zonoiddepthiscontinuousinsidetheconvexhullofthesample,andthusoptimizationw.r.t.
(2)canbeuseddirectly. Byconstructionzonoiddepthcanberepresentedasaproblemoflinear
programming, see Mosler (2002) for details. To account for the missingness, a modification is
necessary, which consists in removing constraints corresponding to the point’s missing coordi-
nates:
min γ s. t. X λ = x , (5)
obs(x) obs
λ(cid:48)1 = 1,
n
γ1 −λ ≥ 0 ,
n n
λ ≥ 0 .
n
Here X the completed n × |obs(x)| data matrix containing columns corresponding only
obs(x)
tononmissingcoordinatesofxand1 (respectively0 )avectorofones(respectivelyzeros)of
n n
lengthn. Imputationisfinallyperformedwiththeλ-weightedaverage
y = X(cid:48) λ.
miss miss(x)
Investigating λ = {λ ,...,λ }, the weights applied to the imputed point, may give additional
1 n
insights concerning its positioning w.r.t. the sample. Thus, from (5) it is clear that in case it
(cid:12)(cid:8) (cid:9)(cid:12)
is solvable, the number of nonnegative weights (cid:12) i : λi > 0, i ∈ {1,...,n} (cid:12) = m + 1,
1 ≤ m+1 ≤ n,whereusuallymofthemareequaltothesolutionγ∗,andoneis≥ 0and≤ γ∗.
Thismeansthaty canbeseenastheaverageofmpointsofthesampleveryslightlyshifted
miss
10