Robust principal components for irregularly spaced longitudinal data
Ricardo A. Maronna





Robust principal components for irregularly spaced longitudinal data
    Ricardo A. Maronna
University of La Plata and University of Buenos Aires (rmaronna@retina.ar)
    December 30, 2023
=================================================================================================




Consider longitudinal data x_ij, with i=1,...,n and j=1,...,p_i,
where x_ij is the j-th observation of the random function X_i(
.) observed at time t_j. The goal of this paper is to develop a
parsimonious representation of the data by a linear combination of a set of
q smooth functions H_k(  .) (k=1,..,q) in the sense that
x_ij≈μ_j+∑_k=1^qβ_kiH_k(  t_j)  ,
such that it fulfills three goals: it is resistant to atypical X_i's
(“case contamination”), it is resistant to
isolated gross errors at some t_ij (”cell
contamination”), and it can be applied when some of the
x_ij are missing (”irregularly spaced —or
'incomplete'– data”).

Two approaches will be proposed for this problem. One deals with the three
goals stated above, and is based on ideas similar to MM-estimation (Yohai
1987). The other is a simple and fast estimator which can be applied to
complete data with case- and cellwise contamination, and is based on applying
a standard robust principal components estimate and smoothing the principal
directions. Experiments with real and simulated data suggest that with
complete data the simple estimator outperforms its competitors, while the MM
estimator is competitive for incomplete data.

Keywords: Principal components, MM-estimator, longitudinal .data, B-splines,
incomplete data.





§ INTRODUCTION


Consider longitudinal data x_ij, with i=1,...,n and j=1,...,p, where
x_ij is the j-th observation of the random function X_i(
.) observed at time t_j. The goal of this paper is to develop a
parsimonious representation of the data by a linear combination of a set of
q smooth functions H_k(  .) (k=1,..,q) in the sense that

    x_ij≈μ_j+∑_k=1^qβ_ikH_k(  t_j)  ,

such that it fulfills three goals: it is resistant to atypical X_i's
(“case contamination”), it is resistant to
isolated gross errors at some t_j (”cell
contamination”), and it can be applied when some of the
x_ij are missing (”irregularly spaced —or
'incomplete'– data”).

Among the abundant literature on this subject, Bali et al. (2011) propose an
approach based on projection pursuit; Boente et al. (2015) and Cevallos (2016)
transform the functional data to lower dimensional data through a spline basis
representation, to which they apply robust principal components based on the
minimization of a robust scale; and Lee et al. (2013) propose a penalized
M-estimator. These approaches hold only for complete data with casewise
contamination. James et al. (2001) and Yao et al (2005) propose estimators for
incomplete data. They require computing sample means and covariances and are
therefore not robust.

Two approaches will be proposed for this problem. One deals with the three
goals stated above, and is based on ideas similar to MM-estimation (Yohai
1987); the other is a very simple and fast estimator which can be applied to
complete data with case- and cellwise contamination. It is based on applying a
standard robust principal components estimate and smoothing the principal
directions, and will be called the “Naive” estimator.

The contents of the paper are as follows. Sections <ref> and
<ref> present the MM- and the Naive estimators. Sections
<ref> and <ref> compare their performances for
complete data to those of former proposals through simulations and a real data
example. Sections <ref> and <ref> compare
the performance of MM for incomplete data to that of a former proposal through
simulations and a real data example. Finally Sections <ref> and
<ref> give the details of the computation of the MM- and Naive estimates.



§ THE “MM” ESTIMATOR



Let B_l(  t) (l=1,...,m) be a basis of B-splines defined on
an interval containing t_1,...,t_p. Call q the desired number of components.

For given α=[α_kl], β=[β_ik] and μ=[μ_j], with
i=1,..,n, j=1,...,p, l=1,..,m and k=1,...,q define

    x_ij(  α.β,μ)  =μ_j
    +∑_k=1^qβ_ikH_k(  t_j)

with

    H_k(  t)  =∑_l=1^mα_klB_l(  t)  .



Then the proposed estimator is given by

    (  a,β,μ)  =min_α,β,μ∑
    _i=1^n∑_j∈ J_iσ_j^2ρ(  x_ij-x_ij(  α,β,μ)
    /σ_j)  .

where σ_j are previously computed local scales and
J_i={j:x_ij is non-missing}. If ρ(  t)  =t^2 and
the data are complete, the estimator coincides with the classical principal
components (PC) estimator. The matrix of principal directions
(“eigenvectors”) is

    𝐄=𝐁α^'∈ R^p× q,

where 𝐁 is the matrix with elements b_jl=B_l(
t_j)  .

This estimator will henceforth be called “MM-estimator”. It is computed iteratively, starting from a
deterministic initial estimator. Since the whole procedure is complex, the
details are postponed to Section <ref>.

A robust measure of “unexplained variance”
 for q components is defined as

    V_q=∑_i=1^n∑_j∈ J_iσ_j^2ρ(
    x_ij-x_ij(  α
    ,β,μ)  /σ_j)  ,

and the “proportion of explained variance”
 is given by

    u_q=1-V_q/V_0,

where V_0 is obtained from (<ref>) with α=0 and
β=0.

The number of spline knots is of the form [p/K], where [.] denotes the
integer part. Choosing K through cross-validation did not yield better
results than using a fixed value. After a series of exploratory simulations it
was decided that K=6 yielded satisfactory results.



§ THE “NAIVE”
 ESTIMATOR


A simple proposal complete data is introduced. It consists of four steps:



  * “Clean” local outliers by applying a
robust smoother to each 𝐱_i and then imputing atypical values.


  * Compute q ordinary robust PC's of the “cleaned” 𝐱_is.


  * Smooth the q eigenvectors using a B-spline basis.


  * Orthogonalize them.


A robust “proportion of unexplained
variance” can be computed as explained in Section
<ref>.

In order to simplify the exposition, the details are given in Section
<ref>.



§ COMPLETE DATA: SIMULATION


Instead of using arbitrary functions for the simulation scenarios, it was
considered more realistic to take them from a real data set. We chose the Low
Resolution Spectroscopy (LRS) data set (Bay 1999), which contains n=531
spectra, on the“red” and “blue” ranges. We deal with the “blue” data, which has p=44.

The first two classical components account respectively for 60% and 33% of
the variability. It was therefore considered that q=2 components was a good choice.

The first four classical eigenvectors 𝐞_k (k=1,..,4) and the
mean vector μ were computed, and a polynomial approximation
was fitted to them, so that they can be obtained for any p. The first two
𝐞_ks and the mean vector are plotted in Figure
<ref>













For given p and q the data are generated as multivariate normal with mean
μ and covariance matrix

    Σ=[  σ_ij]  =∑_k=1^q
    π_k𝐞_k𝐞_k^'+π_0𝐈,

where now μ and 𝐞_k are the original vectors
“scaled” to the interval [1,p] rather
than the original [1,43], and 𝐈 is the p-dimensional identity
matrix. The last term induces noise. The simulations were run for q=2 and 3
and under different configurations of the π^'s.



 §.§ Contamination


Two types of contamination are considered. In case-wise contamination,
a proportion ε_case of cases are replaced by
𝐱_0=K√(λ_1)𝐜,where 𝐜 is
orthogonal to 𝐞_k, k=1,..,q, and λ_1 is the first
eigenvalue of Σ. The factor √(λ_1) simplifies
the choice of the range of K. Two choices for 𝐜 were employed:
one was simply 𝐞_q+1, the (q+1)-th eigenvector of
Σ, and the other was a random direction orthogonal to
𝐞_k, k=1,..,q. Since the qualitative results yielded by both
options were similar, only those corresponding to the first one are shown.

The outlier size K ranges between 0.10 and 3, in order to find the
“most malign” configurations for each estimator.

Figure <ref> shows for p=50 and q=2 three “typical” cases and a case outlier with K=1.2.














In cell-wise contamination, each cell is contaminated at random with
probability ε_cell=1/p. This means that each case has
in average one outlier, but about 63% of the cases have some outliers.
In contaminated cells, x_ij is replaced by x_ij+Kσ_jj. Here K
ranges between 1 and 7.



 §.§ Evaluation


Two issues have to be considered to evaluate an estimator's performance at a
given simulation scenario. The first is the precision with which the PCs
reconstruct the data, which means measuring the differences between the fitted
and the observed values. It was observed that even for “good” data, the distribution of the prediction errors
x_ij-x_ij may have heavy tails. For this reason the fit of
each estimator was evaluated through its prediction Mean Absolute Error (MAE),
which was considered as less sensitive than the Mean Squared Error to some
atypical values.

For a given estimator call x_ij the fitted values. For
case-wise contamination call I the set of non-contaminated rows. Then
define for a data set [x_ij]

    MAE=ave{  |x_ij-x_ij|:i∈ I, j=1,..,p}
    .

For cell-wise contamination the average is taken over all cells:

    MAE=ave{  |x_ij-x_ij
    |:i=1,..,n, j=1,..,p}  ,

the justification being that the estimator should be able to “guess” the true value at the contaminated cell because of
the assumed underlying smoothness of the data.

Finally, for each estimator the average MAE over all Monte Carlo replications
is computed.

The second issue is the estimation of the “true” PCs. To this end, the angle α between the
estimated subspace and the subspace generated by {𝐞_k,k=1,..,q}
was computed. Since it was thought that sin(  α) would be
easier to interpret (because it ranges between 0 and one), the averages of
sin(  α) are reported.



§ SIMULATION RESULTS


The estimators involved in the simulation are: Classic, “Naive”, MM, the S-estimator of Boente and Salibian-Barrera
(2015) using a B-spline basis, and the LTS estimator of Cevallos (2016). The R
code for the S and LTS estimators was kindly supplied by the authors.
Unfortunately it was not possible to obtain the code for the estimator of Lee
et al. (2013).

The dimension was set at p=50. The number of components was chosen as q= 2
and 3. Different configurations for the π_js in (<ref>) were
employed. The results given here are for q=2 and π_1=0.6, π
_2=0.3, π_3=0.05. Since the qualitative results yielded by q=3 and
by other configurations of the π_js are similar, they are not reported
here. The number of replications is 200. It would probably be more realistic
to fix a predetermined proportion of “explained
variance” (say, 90%) and let each estimator choose q
accordingly. While this can be done with the Classic, Naive and MM estimators,
it cannot be done with the S and LTS estimators.

“Efficiency” will be loosely defined as
“similitude with the Classical estimator when ε
_case=ε_cell=0”.

Table <ref> shows the maximum averages (over K) of the MAEs for p=50
and two components.













Table <ref> shows the maximum mean sin(  α) for
the same scenarios.















 §.§ Discussion


Comparing the maximum MAEs it can be concluded that



  * The Naive, S and MM estimators appear as reasonably efficient


  * For case-wise contamination, the Naive and S estimators appear as the
best, the latter being slightly better; LTS and MM have similar acceptable
performances for ε_case=0.1; when ε
_case=0.2, MM and especially LTS have poorer performances


  * For cell-wise contamination, the Naive and MM estimators are practically
unaffected by outliers, while the other estimators are clearly affected by
them.

The results of Table <ref> for ε_cell=0
parallel those of Table <ref>. However, comparing the results for
ε_cell=0.02 with those for ε_cell
=ε_case=0 it is seen that S and LTS are almost
unaffected by cell outliers!. The seeming contradiction between these results
and those of Table <ref> is explained by the fact that in
(<ref>), even if the H_ks are correctly estimated, the βs
must be robustly estimated too, and it is here that these two estimators fail.

Examination of intermediate results reveals that the
“weak spot” of MM lies at the starting
values, which require the estimation of a covariance matrix, which has to be
deterministic, fast, resistant to cell-wise outliers and computable for
incomplete data. The price for so many requirements is a decrease in
robustness for large ε_case. See Section
<ref> for details.



§ COMPLETE DATA: A REAL EXAMPLE


The chosen data set is the “red” range of
the LRS data (Bay 1999), which has n=531 spectra with p=49.

The estimated proportions of explained variance with two components are



c]rrrClassic     Naive     MM

0.91     0.85     0.88







It was therefore considered reasonable to choose q=2 components.

Since our goal is assessing how well do the fitted values approximate the data
(rather than outlier detection), estimates will be evaluated by their MAEs.
For a given estimate let the MAE of row i=1,..,n be

    m_i=ave{| x_ij-x_ij|
    :j=1,...,p}  ,

        

The overall MAEs (averages of the m_i) are

    c]ccccc
    Class.     Naive     MM     S     LTS
    1r184   1r160   1r162   1r160   1r384
       .



Figure <ref> shows the quantiles (in log scale) of the MAEs
corresponding to the different estimates.














The MM, Naive and S estimate give similar results; the classical estimate
gives slightly poorer results, and LTS shows a poor fit. The quantiles of the
classical estimate are larger than those of MM, Naive and S, except for the
top 4%. This fact can be attributed to outlying cases. Figure <ref>
shows three “typical” cases corresponding
to the quartiles of the Naive estimate's m_i (dashed lines) and to the
0.96-quantile (solid line). The latter appears as totally atypical. One may
thus interpret that the classical estimate attempts to fit all cases,
including the atypical ones, at the expense of the typical ones, while the
robust estimates downweight the latter.














For a closer look at the behavior of LTS, Figure <ref> shows the
fit to case 338, a typical one. It is seen that the Naive, S and MM estimators
give good fits and the Classical performs a bit worse, while LTS gives wrong
values at the extremes. This behavior occurs in most cases of this data set.















§ COMPUTING TIMES


The next table shows the average computing times in seconds for n=4p and
number of components q, on a PC with a 3.60 GH processor with 4 cores. Here
“S” denotes the parallelized version of the
S estimator's code.













It is seen that the Naive estimator is by far the fastest, and S is
the slowest; MM is faster than LTS for p≤100, but slower otherwise. The
computing times of S and LTS are unaffected by q; the computing time of MM
increases with q (because components are computed one at a time)

As a final conclusion: considering efficiency, robustness and speed, it may be
concluded that for complete data, the Naive estimator is to be recommended.



§ INCOMPLETE DATA: A REAL EXAMPLE


We begin our treatment of incomplete data by introducing a new data set, which
belongs to the Multicenter AIDS Cohort Study (MACS) and is available at
https://statepi.jhsph.edu/macs/pdt.html. It contains n=305 cases, with
p=48. The number of observations per case varies between 1 and 9. There is a
total of 2001 measurements, which yields a “decimation
rate”

    d=#(measurements)/np=0.137.



The MM-estimator was compared to the procedure “Principal
Analysis by Conditional Estimation” (PACE) (Yao et al 2005),
which is implemented in the package  available at the
authors' web site.

The proportions of explained variance for one and two components were 0.92 and
0.98 for PACE, and 0.77 and 0.81 for MM. It was decided to take q=2.

The overall MAEs of PACE and MM are respectively 93.2 and 87.5. Figure
<ref> compares the quantiles of the MAEs m_i (<ref>)
from both estimates. It is seen that their performances are similar, with MM
doing slightly better.













The proportions of explained variance from MM are much lower than those from
PACE. The reason is that they employ different measures of variability.
Outliers in the first principal axes inflate the ordinary variance but
contribute little to a robust variance. A similar phenomenon is observed in
(Maronna et al, 2006, ps. 213-14).



§ INCOMPLETE DATA: SIMULATION: 


Two scenarios were chosen for the simulations. One is the same employed in
Section <ref>. The other is based on the MACS data set; the mean
vector μ and the first principal directions 𝐞_k of
the data were estimated with PACE, and then employed in the same way as with
the LRS data. They will be called “LRS
scenario” and “MACS
scenario”.

For each scenario a complete data set was generated, and then each cell was
kept at random with probability d, the “decimation
rate” as in (<ref>), and the others set as
“missing”. The values of d were 0.25 and
0.5. The outlier size K ranged between 1 and 6. The values of π_1,
π_2 and π_3 are the same as in Section <ref>. The
dimension was set as p=50, and the sample size n took on the values 100
and 200. Since both values yielded similar results, only those corresponding
to the first one are shown.

Table <ref> shows the results for the MACS scenario with q=2. It
is seen that MM is practically unaffected by the contamination. The MAEs of
PACE are heavily affected by both case- and cell-contamination; the angle is
also heavily affected by case contamination, but not so much by cell contamination.











Figure <ref> shows the MAEs and mean angles as a function of K
for d=0.5. ε_case=0.1 and ε_cell
=0, and Figure <ref> does the same for ε
_case=0 and ε_cell=0.05. The values for
K=0 are those corresponding ε_case=ε
_cell=0. It is seen that in all cases MM remains practically
constant with K; when ε_case=0.1 PACE is clearly
affected when K≥2; and when ε_cell=0.05 its MAEs
are also affected but the angle shows more resistance.



























A surprising feature of Table <ref> is that for ε
_case=ε_cell=0 MM is slightly better than
PACE. Changing the random generator seed yielded similar results.

Table <ref> shows the results for the LRS scenario.











It is seen that when d=0.5 MM is unaffected by contamination and again
outperforms PACE in all cases. When d=0.25 MM has a lower MAE for
ε_case>0. When ε_case=0. PACE
has lower MAE and mean Angle than MM when ε_cell=0, and
a lower maximum mean Angle when ε_cell>0. Since PACE
makes a clever use of the information when data is sparse, this result is more
on line with what would be expected.

It can be concluded that MM is resistant to both case- and cell-contamination,
that its efficiency with respect to PACE is high for decimation rate d=0.5,
while it may depend on the data structure for d=0.25.



§ COMPUTING ALGORITHM OF THE MM-ESTIMATOR


The computation proceeds one component at a time. Define for i=1,..,n and
j=1,...,p

    J_i={j: x_ij is non-missing},  I_j={i: x_ij
     is non-missing}.



At the beginning (“zero components”): apply
robust nonparametric regression to obtain robust and smooth local location and
scale values μ_0j and σ_0j, j=1,...,p as
follows. Let S(  .,.) be a robust smoother. Let m_j be a
location M-estimator of {x_ij, i∈ I_j}; then the set {μ_0j, j=1,..,p} is obtained by applying S to (  τ_j
,m_j,j=1,...,p)  . Let s_j be a τ-scale of {x_ij
-μ_0j,i∈ I_j}; then the set {σ
_0j,j=1,,,,p} is obtained by applying S to (  τ_j
,s_j,j=1,..,p)  . The chosen smoother was the robust version of Loess
(Cleveland 1979) with a span of 0.3.

Let y_ij^(  0)  =x_ij-μ_0j. Compute the
“unexplained variance”

    V_0=1/N∑_i=1^n∑_j∈ I_j^σ_0
    _j^2ρ(  y_ij^(  0)  /σ
    _0_j)    with  N=∑_j=1^pcard(
    I_j)  .



For component 1 use the y_ij^(  0) as input and compute


    (  a^(  1)  ,β^(  1)  ,μ^(  1)
    )  =min_α,β,μ∑_i=1^n
    ∑_j∈ I_jσ_0j^2ρ(  y_ij^(
    0)  -y_ij^(  0)  (  α
    ,β,μ)  /σ_0j)  .



The minimum is computed iteratively, starting from a deterministic initial
estimator to be described in Section <ref>.

Compute the residuals y_ij^(  1)  =y_ij^(  0)
-y_ij^(  0)  (  α
,β,μ)  . Apply a smoother to compute local residual
scales σ_1_j and the “unexplained
variance” with one component:

    V_1=1/N∑_i=1^n∑_j∈ I_jσ_1_j
    ^2ρ(  y_ij^(  1)  /σ_1_j
    )  .

For component k we have

    (  a^(  k)  ,β^(  k)  ,μ^(  k)
    )  =min_α,β,μ∑_i=1^n
    ∑_j∈ I_jσ_k-1,j^2ρ(  y_ij^(
    k-1)  -y_ij^(  k-1)  (
    α,β,μ)  /σ_k-1,j)  .



Each component is orthogonalized with respect to the former ones. The
procedure stops either at a fixed number of components or when the proportion
of explained variance (<ref>) is larger than a given value (e.g. 0.90).



 §.§ The iterative algorithm


Computing each component requires an iterative algorithm and starting values.
The algorithm is essentially one of “alternating
regressions”.

Recall that at each step α∈R^m and β∈R^n are one-dimensional. Put as usual ψ=ρ^' and
W(s)=ψ(  s)  /s. Put for brevity h(  t)
=∑_l=1^mα_lB_l(  t) where B_i are the
elements of the spline basis.

Differentiating the criterion in (<ref>) yields a set of estimating
equations that can be written in fixed-point form, yielding a
“weighted alternating regressions” scheme.
To simplify the notation the superscript (  k-1) will be
dropped from y_ij^(  k-1) and y_ij^(
k-1)  . Put

    w_ij=w_ij(  α,β,μ)  =W(
    y_ij-y_ij(  α,β,μ)
    /σ_k-1,j)  .

Then μ_j and β_i can be expressed as weighted residual means and
weighted univariate least squares regressions, respectively:

    μ_j=1/∑_i∈ I_jw_ij∑_i∈ I_jw_ij(
    y_ij-β_ih(  t_j)  )  ,


    β_i=∑_j∈ J_iw_ijh(  t_j)  (
    y_ij-μ_j)/∑_j∈ J_iw_ijh(  t_j)  ^2

and α is the solution of

    ∑_i=1^n∑_j∈ J_iw_ij(  y_ij-μ_j)  β
    _i𝐛(  t_j)  =∑_i=1^n∑_j∈ J_i
    w_ijβ_i^2𝐛(  t_j)  𝐛(
    t_j)  ^'α

with 𝐛(  t)  =(  B_1(  t)
,...,B_l(  t)  )'.

At each iteration the w_ij are updated. It can be shown that the criterion
descends at each iteration.



 §.§ The initial values


For each component, initial values for α and
β are needed. They should be deterministic, since
subsampling would make the procedure impractically slow.



  §.§.§ The initial α


For k,l∈{  1,...,p} call N_kl the number of cases that
have values in both t_k and t_l:

    N_kl=#(  I_k∩ I_l)  .



In longitudinal studies, many N_kl may null or very small.

Compute a (possibly incomplete) p× p matrix Σ
=[σ_kl] of pairwise robust covariances of (
y_ik,y_il:i∈ I_k∩ I_l) with the Gnanadesikan-Kettenring
procedure:

    Cov(  X,Y)  =1/4(  S(  X+Y)
    ^2-S(  X-Y)  ^2)  ,

where S is a robust dispersion. Preliminary simulations led to the choice of
the Q_n estimator of Rousseeuw and Croux (1993).

More precisely; for N_kl≥3 compute σ_kl as above. If
min_klN_kl is “large enough” (here:
≥10) use the resulting Σ.

Otherwise apply a two-dimensional smoother to improve Σ and to
fill in the missing values. The bivariate Loess was employed for this purpose.

Then compute the first eigenvector 𝐞 of Σ (note
that Σ is not guaranteed to be positive definite and hence
further principal components may be unreliable). Given 𝐞,
smooth it using the spline basis. This α follows from
(<ref>).



  §.§.§ The initial β


For i=1,...,n the initial β_i is a robust univariate regression of
y_ij on h(  t_j) (j∈ J_i), namely the L_1
regression, which is fast and reliable.

Note that only cellwise outliers matter at this step.



 §.§ The final adjustment


Note that the former steps yield only an approximate solution of
(<ref>) since the components are computed one at a time. In order to
improve the approximation a natural procedure is as follows. After computing
q components we have a p× q-matrix 𝐔 of principal
directions with elements

    u_jk=∑_l=1^mB_l(  t_j)  α_kl,

an n× p-matrix of weights 𝐖. and a location vector
μ. Then a natural improvement is –keeping 𝐔 and
μ fixed– to recompute the βs by means of univariate
weighted regressions with weights w_ij. Let β_i.
=[β_ik,k=1,..,p], and set

    β_i=min_β∈R^q∑_j∈ J_i
    w_ij(  x_ij-μ_j-𝐔β)  ^2.

The effect of this step in the case of complete data is negligible, but it
does improve the estimator's behavior for incomplete data.

However, it was found out that the improvement is not good enough when the
data are very sparse. For this reason another approach was used, namely, to
compute β_i as a regression M-estimate. Let 𝐳
_i=(x_ij-μ_j:j∈ J_i) and 𝐕=[v_jk] with
v_jk=u_jk for j∈ J_i. Then β_i is a bisquare
regression estimate of 𝐳_i on 𝐕, with tuning constant
equal to 4, using L_1 as a starting estimate. Note that here only cell
outliers matter, and therefore L_1 yields reliable starting values. The
estimator resulting from this step does not necessarily coincide with
(<ref>), but simulations show that it is much better than the
“natural” adjustment described above when
the data are very sparse.



§ THE “NAIVE” ESTIMATOR:
DETAILS


In step 1 of Section <ref>, compute for each 𝐱_i
robust local location and scatter estimates μ_i
,σ_i. The “cleaned”
 values are

    x_ij=μ_i+σ_iψ(
    x_ij-μ_i/σ_i)  ,

where ψ is the bisquare ψ-function with tuning constant equal to 4.

The ordinary robust PCs of step 2 are computed using the cleaned data
x_ij with the S-M estimator of (Maronna 2005). Call
{𝐱_i^(  q)  } the fit for q components
and put r_i^(  q)  = ‖𝐱_i
-𝐱_i^(  q)  ‖ , i=1,...,n. Then
the estimator minimizes S(  r_i^(  q)  
,  i=1,..,n) where S is the bisquare M-scale. The
“proportion of unexplained variance” is

    S(  r_i^(  q)  ,  i=1,..,n)  /S(
    r_i^(  0)  ,  i=1,..,n)  .



The number of knots in step 3 is chosen through generalized cross-validation.



§ REFERENCES


Bali, J.L., Boente, G., Tyler, D.E. and Wang, J-L. (2011). Robust functional
principal components: a projection-pursuit approach. The Annals of
Statistics, 39, 2852–2882.

Bay, S. D. (1999), The UCI KDD Archive [http://kdd.ics.uci.edu], University of
California, Irvine, Dept. of Information and Computer Science.

Boente, G. and Salibian-Barrera, M. (2015). S-Estimators for Functional
Principal Component Analysis. Journal of the American Statistical
Association, 110, 1100-1111.

Cevallos Valdiviezo, H. (2016). On Methods for Prediction Based on Complex
Data with Missing Values and Robust Principal Component Analysis, PhD thesis,
Ghent University (supervisors Van Aelst S. and Van den Poel, D.).

Cleveland, W.S. (1979). Robust Locally Weighted Regression and Smoothing
Scatterplots. Journal of the American Statistical Association,
74, 829-836.

James, G., Hastie, T.G., and Sugar, C.A. (2001). Principal Component Models
for Sparse Functional Data. Biometrika, 87, 587-602.

Lee, S., Shin, H. and Billor, N. (2013). M-type smoothing spline estimators
for principal functions. Computational Statistics and Data Analysis,
66, 89-100.

Maronna, R. (2005). Principal components and orthogonal regression based on
robust scales. Technometrics, 47




264-273.

Maronna, R.A., Martin, R.D. and Yohai, V.J. (2006). Robust
Statistics: Theory and Methods. John Wiley and Sons, Chichester.

Rousseeuw, P.J. and Croux, C. (1993). Alternatives to the Median Absolute
Deviation. Journal of the American Statistical Association,
88, 1273–1283.

Yao, F., Müller, H-G. and Wang, J-L. (2005). Functional Data Analysis for
Sparse Longitudinal Data. Journal of the American Statistical
Association, 100, 577-590.

Yohai, V.J. (1987). High Breakdown-Point and High Efficiency Robust Estimates
for Regression. The Annals of Statistics, 15, 642-656.


