# Robust Multivariate Functional Control Charts

Christian Capezza  Fabio Centofanti Corresponding author. e-mail:
<EMAIL_ADDRESS>Antonio Lepore  Biagio Palumbo

###### Abstract

Profile monitoring assesses the stability over time of one or multiple
functional quality characteristics to quickly detect special causes of
variation that act on a process. In modern Industry 4.0 applications, a huge
amount of data is acquired during manufacturing processes that are often
contaminated with anomalous observations in the form of both casewise and
cellwise outliers. Because anomalous observations can seriously affect the
monitoring performance, profile monitoring methods that are able to
successfully deal with outliers are needed. To this aim, we propose a new
framework, referred to as robust multivariate functional control charts
(RoMFCC), that is able to monitor multivariate functional data while being
robust to both functional casewise and cellwise outliers. The RoMFCC relies on
four main elements: (I) a univariate filter to identify functional cellwise
outliers to be replaced by missing components; (II) a robust functional data
imputation method of missing values; (III) a casewise robust dimensionality
reduction; (IV) a monitoring strategy for the multivariate functional quality
characteristic. An extensive Monte Carlo simulation study is performed to
quantify the monitoring performance of the RoMFCC by comparing it with some
competing methods already available in the literature. Finally, a motivating
real-case study is presented where the proposed framework is used to monitor a
resistance spot welding process in the automotive industry.

Keywords: Functional Data Analysis, Profile Monitoring, Statistical Process
Control, Robust Estimation, Casewise and Cellwise Outliers

##  1 Introduction

In modern industrial applications, data acquisition systems allow to collect
massive amounts of data with high frequency. Several examples may be found in
the current Industry 4.0 framework, which is reshaping the variety of signals
and measurements that can be gathered during manufacturing processes. The
focus in many of these applications is statistical process monitoring (SPM),
whose main aim is to quickly detect unusual conditions in a process when
special causes of variation act on it, i.e., the process is out-of-control
(OC). On the contrary, when only common causes are present, the process is
said to be in-control (IC).

In this context, the experimental measurements of the quality characteristic
of interest are often characterized by complex and high dimensional formats
that can be well represented as functional data or profiles (Ramsay and
Silverman, , 2005; Kokoszka and Reimherr, , 2017). The simplest approach for
monitoring one or multiple functional variables is based on the extraction of
scalar features from each profile, e.g., the mean, followed by the application
of classical SPM techniques for multivariate data (Montgomery, , 2012).
However, feature extraction is known to be problem-specific, arbitrary, and
risks compressing useful information. Thus, there is a growing interest in
profile monitoring (Noorossana et al., , 2011), whose aim is to monitor a
process when the quality characteristic is best characterized by one or
multiple profiles. Some recent examples of profile monitoring applications can
be found in Menafoglio et al.,  (2018); Capezza et al.,  (2020); Capezza et
al., 2021a ; Capezza et al., 2021b ; Centofanti et al., 2021b .

The main tools for SPM are control charts that are implemented in two phases.
In Phase I, historical process data are used to set control chart limits to be
used in Phase II, i.e., the actual monitoring phase, where observations
falling outside the control limits are signaled as OC. In classical SPM
applications the historical Phase I data are assumed to come from an IC
process. However, this assumption is not always valid. As an example, let
consider the motivating real-world application, detailed in Section 4, that
concerns the SPM of a resistance spot welding (RSW) process in the assembly of
automobiles. RSW is the most common technique employed in joining metal sheets
during body-in-white assembly of automobiles (Zhou and Cai, , 2014), mainly
because of its adaptability for mass production (Martín et al., , 2014). Among
on-line measurements of RSW process parameters, the so-called dynamic
resistance curve (DRC) is recognized as the full technological signature of
the metallurgical development of a spot weld (Dickinson et al., , 1980;
Capezza et al., 2021b, ) and, thus, can be used to characterize the quality of
a finished sub-assembly. Figure 1 shows 100 DRCs corresponding to 10 different
spot welds, measured in $m\Omega$, that are acquired during the RSW process
from the real-case study presented in Section 4. Several outliers are clearly
visible that should be taken into account to set up an effective SPM strategy.

Figure 1: Sample of 100 DRCs, measured in $m\Omega$, that are acquired during
the RSW process from the real-case study in Section 4. The different panels
refer to the corresponding different spot welds, denoted with names from V1 to
V10.

Indeed, control charts are very sensitive to the presence of outlying
observations in Phase I that can lead to inﬂated control limits and reduced
power to detect process changes in Phase II. To deal with outliers in the
Phase I sample, SPM methods use two common alternatives, namely the diagnostic
and the robust approaches (Kruger and Xie, , 2012; Hubert et al., , 2015). The
diagnostic approach is based on standard estimates after the removal of sample
units identified as outliers that translates into SPM methods where iterative
re-estimation procedures are considered in Phase I. This approach could be
often safely applied to eliminate the eﬀect of a small number of very extreme
observations. However, it will fail to detect more moderate outliers that are
not always as easy to label correctly. On the contrary, the robust approach
accepts all the data points and tries to ﬁnd a robust estimator which reduces
the impact of outliers on the final results (Maronna et al., , 2019).

Several robust approaches for the SPM of a multivariate scalar quality
characteristic have been proposed in the literature. Alfaro and Ortega,
(2009) show a comparison of robust alternatives to the classical Hotelling’s
control chart. It includes two alternative Hotelling’s $T^{2}$-type control
charts for individual observations, proposed by Vargas,  (2003) and Jensen et
al.,  (2007), which are based on the minimum volume ellipsoide and the minimum
covariance determinant estimators (Rousseeuw, , 1984), respectively. Moreover,
the comparison includes the control chart based on the reweighted minimum
covariance determinant (RMCD) estimators, proposed by Chenouri et al.,
(2009). More recently, Cabana and Lillo,  (2021) propose an alternative robust
Hotelling’s $T^{2}$ procedure using the robust shrinkage reweighted estimator.

Although Kordestani et al.,  (2020) and Moheghi et al.,  (2020) propose robust
estimators to monitoring simple linear profiles, to the best of authors’
knowledge, a robust approach able to successfully capture the functional
nature of a multivariate functional quality characteristic has not been
devised in the SPM literature so far. Beyond the SPM literature, several works
have been proposed to deal with outlying functional observations. Several
methods extend the classical linear combination type estimators (i.e.,
L-estimator) (Maronna et al., , 2019) to the functional setting to robustly
estimate the center of a functional distribution through trimming (Fraiman and
Muniz, , 2001; Cuesta-Albertos and Fraiman, , 2006) and functional data depths
(Cuevas and Fraiman, , 2009; López-Pintado and Romo, , 2011). Sinova et al.,
(2018) introduce the notion of maximum likelihood type estimators (i.e.,
M-estimators) in the functional data setting. More recently, Centofanti et
al., 2021a propose a robust functional ANOVA method that reduces the weights
of outlying functional data on the results of the analysis.

Regarding functional principal component analysis (FPCA), robust approaches
are classified by Boente and Salibián-Barrera,  (2021) in three groups,
depending on the specific property of principal components on which they
focus. Methods in the first group perform the eigenanalysis of a robust
estimator of the scatter operator, as the spherical principal components
method of Locantore et al.,  (1999) and the indirect approach of Sawant et
al.,  (2012). The latter performs a robust PCA method, e.g., ROBPCA (Hubert et
al., , 2005), on the matrix of the basis coefficients corresponding to a basis
expansion representation of the functional data. The second group includes
projection-pursuit approaches (Hyndman and Ullah, , 2007), which sequentially
search for the directions that maximize a robust estimator of the spread of
the data projections. Whereas, the third group is composed of methods that
estimate the principal components spaces by minimizing a robust reconstruction
error measure (Lee et al., , 2013). Finally, it is worth mentioning diagnostic
approaches for functional outlier detection, which have been proposed for both
univariate (Hyndman and Shang, , 2010; Arribas-Gil and Romo, , 2014; Febrero
et al., , 2008) and multivariate functional data (Hubert et al., , 2015; Dai
and Genton, , 2018; Alemán-Gómez et al., , 2022).

In presence of many functional variables, the lack of robust approaches that
deal with outliers is exacerbated by the curse of dimensionality. Traditional
multivariate robust estimators assume a casewise contamination model for the
data, which consists in a mixture of two distributions, where the majority of
the cases is free of contamination and the minority mixture component
describes an unspecified outlier generating distribution. Alqallaf et al.,
(2009) show that these traditional estimators are affected by the problem of
propagation of outliers. In situations where only a small number of cases are
contaminated, the traditional robust approaches work well. However, under an
independent contamination model such as cellwise outliers (i.e., contamination
in each variable is independent from the other variables), when the
dimensionality of the data is high, the fraction of perfectly observed cases
can be rather small and the traditional robust estimators may fail. Moreover,
Agostinelli et al.,  (2015) point out that both types of data contamination,
casewise and cellwise, may occur together. This problem has been addressed in
the multivariate scalar setting by Agostinelli et al.,  (2015) that propose a
two steps method. In the first step, a univariate filter is used to eliminate
large cellwise outliers, i.e., detection and replacement by missing values,
then, in the second step, a robust estimation, specifically designed to deal
with missing data, is applied to the incomplete data. Leung et al.,  (2016)
notice that the univariate filter does not handle well moderate-size cellwise
outliers, therefore they introduce for the first step a consistent bivariate
filter to be used in combination with the univariate filter. Rousseeuw and
Bossche,  (2018) propose a method for detecting deviating data cells that
takes the correlations between the variables into account, whereas, Tarr et
al.,  (2016) devise a method for robust estimation of precision matrices under
cellwise contamination. Other methods that consider cellwise outliers have
been developed for regression and classification (Filzmoser et al., , 2020;
Aerts and Wilms, , 2017).

To deal with multivariate functional outliers, in this paper we propose a new
framework, referred to as robust multivariate functional control chart
(RoMFCC), for SPM of multivariate functional data that is robust to both
functional casewise and cellwise outliers. The latter corresponds to a
contamination model where outliers arise in each variable independently from
the other functional variables. Specifically, to deal with functional cellwise
outliers, the proposed framework considers an extension of the filtering
approach proposed by Agostinelli et al.,  (2015) to univariate functional data
and an imputation method inspired by the robust imputation technique of
Branden and Verboven,  (2009). Moreover, it also considers a robust
multivariate functional principal component analysis (RoMFPCA) based on the
ROBPCA method (Hubert et al., , 2005), and a profile monitoring strategy built
on the Hotelling’s $T^{2}$ and the squared prediction error ($SPE$) control
charts (Noorossana et al., , 2011; Grasso et al., , 2016; Centofanti et al.,
2021b, ; Capezza et al., , 2020; Capezza et al., 2021a, ; Capezza et al., ,
2022). A Monte Carlo simulation study is performed to quantify the probability
of signal (i.e., detecting an OC observation) of RoMFCC in identifying mean
shifts in the functional variables in presence of both casewise and cellwise
outliers. That is, the proposed RoMFCC is compared with other control charts
already present in the literature. Finally, the practical applicability of the
proposed control chart is illustrated on the motivating real-case study in the
automotive manufacturing industry. In particular, the RoMFCC is shown to
adequately identify a drift in the manufacturing process due to electrode
wear.

The article is structured as follows. Section 2 introduces the proposed RoMFCC
framework. In Section 3, a simulation study is presented where RoMFCC is
compared to other popular competing methods. The real-case study in the
automotive industry is presented in Section 4. Section 5 concludes the
article. Supplementary materials for this article are available online. All
computations and plots have been obtained using the programming language R (R
Core Team, , 2021).

##  2 The Robust Multivariate Functional Control Chart Framework

The proposed RoMFCC is a new general framework for SPM of multivariate
functional data that is able to deal with both functional casewise and
cellwise outliers. It relies on the following four main elements.

  1. (I)

Functional univariate filter, which is used to identify functional cellwise
outliers to be replaced by missing components.

  2. (II)

Robust functional data imputation, where a robust imputation method is applied
to the incomplete data to replace missing values.

  3. (III)

Casewise robust dimensionality reduction, which reduces the infinite
dimensionality of the multivariate functional data by being robust towards
casewise outliers.

  4. (IV)

Monitoring strategy, to appropriately monitor multivariate functional data.

In what follows, we describe a specific implementation of the RoMFCC framework
where (I) an extension of the filtering proposed by Agostinelli et al.,
(2015), referred to as functional univariate filter (FUF), is considered; (II)
a robust functional data imputation method referred to as RoFDI and based on
the robust imputation technique of Branden and Verboven,  (2009) is used;
(III) the RoMFPCA is considered as the casewise robust dimensionality
reduction method. Finally, (IV) the multivariate functional data are monitored
through the profile monitoring approach based on the simultaneous application
of the Hotelling’s $T^{2}$ and the squared prediction error ($SPE$) control
charts. For ease of presentation, the RoMFPCA is presented in Section 2.1.
Then, Section 2.2, Section 2.3, and, Section 2.4 describe the FUF, the RoFDI
method, and the monitoring strategy, respectively. Section 2.5 details the
Phase I and Phase II of the proposed implementation of the RoMFCC framework
where the elements (I-IV) are put together.

###  2.1 Robust Multivariate Functional Principal Component Analysis

Let $\bm{X}=\left(X_{1},\dots,X_{p}\right)^{T}$ a random vector with
realization in the Hilbert space $\mathbb{H}$ of $p$-dimensional vectors of
functions defined on the compact set $\mathcal{T}\in\mathbb{R}$ with
realizations in $L^{2}(\mathcal{T})$, i.e., the Hilbert spaces of square
integrable functions defined on $\mathcal{T}$. Accordingly, the inner product
of two functions $f$ and $g$ in $L^{2}\left(\mathcal{T}\right)$ is $\langle
f,g\rangle=\int_{\mathcal{T}}f\left(t\right)g\left(t\right)dt$, and the norm
is $\lVert\cdot\rVert=\sqrt{\langle\cdot,\cdot\rangle}$. The inner product of
two function vectors $\mathbf{f}=\left(f_{1},\dots,f_{p}\right)^{T}$ and
$\mathbf{g}=\left(g_{1},\dots,g_{p}\right)^{T}$ in $\mathbb{H}$ is
$\langle\mathbf{f},\mathbf{g}\rangle_{\mathbb{H}}=\sum_{j=1}^{p}\langle
f_{j},g_{j}\rangle$ and the norm is
$\lVert\cdot\rVert_{\mathbb{H}}=\sqrt{\langle\cdot,\cdot\rangle_{\mathbb{H}}}$.

We assume that $\bm{X}$ has mean
$\bm{\mu}=\left(\mu_{1},\dots,\mu_{p}\right)^{T}$,
$\mu_{i}(t)=\operatorname{E}(X_{i}(t))$, $t\in\mathcal{T}$ and covariance
$\bm{G}=\\{G_{ij}\\}_{1\leq i,j\leq p}$,
$G_{ij}(s,t)=\operatorname{Cov}(X_{j}(s),X_{j}(t))$, $s,t\in\mathcal{T}$. In
what follows, to take into account for differences in degrees of variability
and units of measurements among $X_{1},\dots,X_{p}$, the transformation
approach of Chiou et al.,  (2014) is considered. Specifically, let consider
the vector of standardized variables
$\bm{Z}=\left(Z_{1},\dots,Z_{p}\right)^{T}$,
$Z_{i}(t)=v_{i}(t)^{-1/2}(X_{i}(t)-\mu_{i}(t))$, with $v_{i}(t)=G_{ii}(t,t)$,
$t\in\mathcal{T}$. Then, from the multivariate Karhunen-Loève’s Theorem (Happ
and Greven, , 2018) follows that

$\bm{Z}(t)=\sum_{l=1}^{\infty}\xi_{l}\bm{\psi}_{l}(t),\quad t\in\mathcal{T},$

where $\xi_{l}=\langle\bm{\psi}_{l},\bm{Z}\rangle_{\mathbb{H}}$ are random
variables, said principal components scores or simply scores such that
$\operatorname{E}\left(\xi_{l}\right)=0$ and
$\operatorname{E}\left(\xi_{l}\xi_{m}\right)=\lambda_{l}\delta_{lm}$, with
$\delta_{lm}$ the Kronecker delta. The elements of the orthonormal set
$\\{\bm{\psi}_{l}\\}$,
$\bm{\psi}_{l}=\left(\psi_{l1},\dots,\psi_{lp}\right)^{T}$, with
$\langle\bm{\psi}_{l},\bm{\psi}_{m}\rangle_{\mathbb{H}}=\delta_{lm}$, are
referred to as principal components, and are the eigenfunctions of the
covariance $\bm{C}$ of $\bm{Z}$ corresponding to the eigenvalues
$\lambda_{1}\geq\lambda_{2}\geq\dots\geq 0$. Following the approach of Ramsay
and Silverman,  (2005), the eigenfunctions and eigenvalues of the covariance
$\bm{C}$ are estimated through a basis function expansion approach.
Specifically, we assume that the functions $Z_{j}$ and a generic eigenfunction
$\bm{\psi}$ of $\bm{C}$ with components $\psi_{j}$, for $j=1,\dots,p$, can be
represented as

$Z_{j}(t)\approx\sum_{k=1}^{K}c_{jk}\phi_{jk}(t),\quad\psi_{j}(t)\approx\sum_{k=1}^{K}b_{jk}\phi_{jk}(t),\quad
t\in\mathcal{T}$ (1)

where $\bm{\phi}_{j}=\left(\phi_{j1},\dots,\phi_{jK}\right)^{T}$,
$\bm{c}_{j}=\left(c_{j1},\dots,c_{jK}\right)^{T}$ and
$\bm{b}_{j}=\left(b_{j1},\dots,b_{jK}\right)^{T}$ are the basis functions and
coefficient vectors for the expansion of $Z_{j}$ and $\psi_{j}$, respectively.
With these assumptions, standard multivariate functional principal component
analysis (Ramsay and Silverman, , 2005; Chiou et al., , 2014) estimates
eigenfunctions and eigenvalues of the covariance $\bm{C}$ by performing
standard multivariate principal component analysis on the random vector
$\bm{W}^{1/2}\bm{c}$, where
$\bm{c}=\left(\bm{c}_{1}^{T},\dots,\bm{c}_{p}^{T}\right)^{T}$ and $\bm{W}$ is
a block-diagonal matrix with diagonal blocks $\bm{W}_{j}$, $j=1,\dots,p$,
whose entries are $w_{k_{1}k_{2}}=\langle\phi_{k_{1}},\phi_{k_{2}}\rangle$,
$k_{1},k_{2}=1,\dots,K$. Then, the eigenvalues of $\bm{C}$ are estimated by
those of the covariance matrix of $\bm{W}^{1/2}\bm{c}$, whereas, the
components $\psi_{j}$ of the generic eigenfunction, with corresponding
eigenvalue $\lambda$, are estimated through Equation (1) with
$\bm{b}_{j}=\bm{W}^{-1/2}\bm{u}_{j}$, where
$\bm{u}=\left(\bm{u}_{1}^{T},\dots,\bm{u}_{p}^{T}\right)^{T}$ is the
eigenvector of the covariance matrix of $\bm{W}^{1/2}\bm{c}$ corresponding to
$\lambda$. However, it is well known that standard multivariate principal
component analysis is not robust to outliers (Maronna et al., , 2019), which
obviously reflects on the functional principal component analysis by probably
providing misleading results. Extending the approach of Sawant et al.,  (2012)
for multivariate functional data, the proposed RoMFPCA applies a robust
principal component analysis alternative to the random vector
$\bm{W}^{1/2}\bm{c}$. Specifically, we consider the ROBPCA approach (Hubert et
al., , 2005), which is a computationally efficient method explicitly conceived
to produce estimates with high breakdown in high dimensional data settings,
which almost always arise in the functional context, and to handle large
percentage of contamination. Thus, given $n$ independent realizations
$\bm{X}_{i}$ of $\bm{X}$, dimensionality reduction is achieved by
approximating $\bm{X}_{i}$ through $\hat{\bm{X}}_{i}$, for $i=1,\dots,n$, as

$\hat{\bm{X}}_{i}(t)=\hat{\bm{\mu}}(t)+\hat{\bm{D}}(t)\sum_{l=1}^{L}\hat{\xi}_{il}\hat{\bm{\psi}}_{l}(t)\quad
t\in\mathcal{T}$ (2)

where $\hat{\bm{D}}$ is a diagonal matrix whose diagonal entries are robust
estimates $\hat{v}_{j}^{1/2}$ of $v_{j}^{1/2}$,
$\hat{\bm{\mu}}=\left(\hat{\mu}_{1},\dots,\hat{\mu}_{p}\right)^{T}$ is a
robust estimate of $\bm{\mu}$, $\hat{\bm{\psi}}_{l}$ are the first $L$
robustly estimated principal components and
$\hat{\xi}_{il}=\langle\hat{\bm{\psi}}_{l},\hat{\bm{Z}}_{i}\rangle_{\mathbb{H}}$
are the estimated scores with robustly estimated variances
$\hat{\lambda}_{l}$. The estimates $\hat{\bm{\psi}}_{l}$ and
$\hat{\lambda}_{l}$ are obtained through the $n$ realizations of $\bm{Z}_{i}$
estimated by using $\hat{\mu}_{j}$ and $\hat{v}_{j}$. The robust estimates
$\hat{\mu}_{j}$ and $\hat{v}_{j}$ are obtained through the scale equivariant
functional $M$-estimator and the functional normalized median absolute
deviation estimator proposed by Centofanti et al., 2021a . As in the
multivariate setting, $L$ is generally chosen such that the retained principal
components $\hat{\bm{\psi}}_{1},\dots,\hat{\bm{\psi}}_{L}$ explain at least a
given percentage of the total variability, which is usually in the range
70-90$\%$, however, more sophisticated methods could be used as well, see
Jolliffe,  (2011) for further details.

###  2.2 Functional Univariate Filter

To extend the filter of Agostinelli et al.,  (2015) and Leung et al.,  (2016)
to univariate functional data, let consider $n$ independent realizations
$X_{i}$ of a random function $X$ with values in $L^{2}(\mathcal{T})$. The
proposed FUF considers the functional distances $D_{i}^{fil}$, $i=1,\dots,n$,
defined as

$D_{i}^{fil}=\sum_{l=1}^{L^{fil}}\frac{(\hat{\xi}_{il}^{fil})^{2}}{\hat{\lambda}_{l}^{fil}},$
(3)

where the estimated scores
$\hat{\xi}_{il}^{fil}=\langle\hat{{\psi}}_{l}^{fil},\hat{{Z}}_{i}\rangle$, the
estimated eigenvalues $\hat{\lambda}_{l}^{fil}$, the estimated principal
components $\hat{{\psi}}_{j}^{fil}$, and the estimated standardized
observations $\hat{{Z}}_{i}$ of ${X}_{i}$ are obtained by applying, with
$p=1$, the RoMFPCA described in Section 2.1 to the sample $X_{i}$,
$i=1,\dots,n$. In this setting, RoMFPCA is used to appropriately represent
distances among $X_{i}$’s and not to perform dimensionality reduction, this
means that $L^{fil}$ should be sufficiently large to capture a large
percentage of the total variability $\delta^{fil}$. Let $G_{n}$ the empirical
distribution of $D_{i}^{fil}$, that is

$G_{n}(x)=\frac{1}{n}\sum_{i=1}^{n}I(D_{i}^{fil}\leq x),\quad x\geq 0,$

where $I$ is the indicator function. Then, functional observations $X_{i}$ are
labeled as cellwise outliers by comparing $G_{n}(x)$ with $G(x)$, $x\geq 0$,
where $G$ is a reference distribution for $D_{i}^{fil}$. Following Leung et
al.,  (2016), we consider the chi-squared distribution with $L^{fil}$ degrees
of freedom, i.e., $G=\chi^{2}_{L^{fil}}$. The proportion of flagged cellwise
outliers is defined by

$d_{n}=\sup_{x\geq\eta}\\{G(x)-G_{n}(x)\\}^{+},$

where $\\{a\\}^{+}$ represents the positive part of $a$, and
$\eta=G^{-1}(\alpha)$ is a large quantile of $G$. Following Agostinelli et
al.,  (2015), in the following we consider $\alpha=0.95$ as the aim is to
detect extreme cellwise outliers, but other choices could be considered as
well. Finally, we flag $\left[nd_{n}\right]$ observations with the largest
functional distances $D_{i}^{fil}$ as functional cellwise outliers (here
$\left[a\right]$ is the largest integer less than or equal to $a$). From the
arguments in Agostinelli et al.,  (2015) and Leung et al.,  (2016) the FUF is
a consistent even when the actual distribution of $D_{i}^{fil}$ is unknown,
that is asymptotically the filter will not wrongly flag a cellwise outlier
provided that the tail of the chosen reference distribution $G$ is heavier (or
equal) than that of the actual unknown distribution.

###  2.3 Robust Functional Data Imputation

Let consider $n$ independent realizations
$\bm{X}_{i}=\left(X_{i1},\dots,X_{ip}\right)^{T}$ of a random vector of
functions $\bm{X}$ as defined in Section 2.1. This section considers the
setting where $\bm{X}_{i}$, $i=1,\dots,n$, may presents missing components,
i.e., at least one of $X_{i1},\dots,X_{ip}$ is missing. Thus, for each
realization $\bm{X}_{i}$ we can identify the missing
$\bm{X}^{m}_{i}=\left(X_{im_{i1}},\dots,X_{im_{is_{i}}}\right)^{T}$ and
observed $\bm{X}^{o}_{i}=\left(X_{io_{i1}},\dots,X_{io_{it_{i}}}\right)^{T}$
components with $t_{i}=p-s_{i}$, and where $\\{m_{ij}\\}$ and $\\{o_{ij}\\}$
are disjoint set of indices, whose union coincides with the set of indeces
$\\{1,\dots,p\\}$, that indicate which components of the realization $i$ are
either observed or missing. Moreover, we assume that a set $S_{c}$ of $c$
realizations free of missing components is available. The proposed RoFDI
method extends to the functional setting the robust imputation approach of
Branden and Verboven,  (2009) that sequentially estimates the missing part of
an incomplete observation such that the imputed observation has minimum
distance from the space generated by the complete realizations. Analogously,
starting from $S_{c}$, we propose that the missing components of the
observation $\bm{X}_{\underline{i}}\notin S_{c}$, which correspond to the
smallest $s_{i}$, are sequentially imputed by minimizing, with respect to
$\bm{X}^{m}_{\underline{i}}$,

$D(\bm{X}_{\underline{i}}^{m})=\sum_{l=1}^{L^{imp}}\frac{(\hat{\xi}_{il}^{imp})^{2}}{\hat{\lambda}_{l}^{imp}},$
(4)

where the estimated scores
$\hat{\xi}_{il}^{imp}=\langle\hat{\bm{\psi}}_{l}^{imp},\hat{\bm{Z}}_{\underline{i}}\rangle_{\mathbb{H}}$,
eigenvalues $\hat{\lambda}_{l}^{imp}$, principal components
$\hat{\bm{\psi}}_{l}^{imp}$ and standardized observations
$\hat{\bm{Z}}_{\underline{i}}$ of $\bm{X}_{\underline{i}}$ are obtained by
applying the RoMFPCA (Section 2.1) to the free of missing data observations in
$S_{c}$. Analogously to Section 2.2, RoMFPCA is used to define the distance of
$\bm{X}_{\underline{i}}$ from the space generated by the free of missing data
observations, thus, $L^{imp}$ should be sufficiently large to capture a large
percentage of the total variability $\delta^{imp}$. Because
$\hat{\bm{Z}}_{\underline{i}}$ is the standardized version of
$\bm{X}_{\underline{i}}$, we can identify the missing
$\hat{\bm{Z}}^{m}_{\underline{i}}$ and observed
$\hat{\bm{Z}}^{o}_{\underline{i}}$ components of
$\hat{\bm{Z}}_{\underline{i}}$. Thus, the minimization problem in Equation (4)
can be equivalently solved with respect to $\hat{\bm{Z}}_{\underline{i}}^{m}$,
and the resulting solution can be unstandardized to obtain the imputed
components of $\bm{X}_{\underline{i}}^{m}$.

Due to the approximations in Equation (1), $\hat{\bm{Z}}_{\underline{i}}$ is
uniquely identified by the coefficient vectors $\bm{c}_{\underline{i}j}$,
$j=1,\dots,p$, related to the basis expansions of its components. Let
$\bm{c}^{m}_{\underline{i}j}$,
$j=m_{\underline{i}1},\dots,m_{\underline{i}s_{\underline{i}}}$ and
$\bm{c}^{o}_{\underline{i}j}$,
$j=o_{\underline{i}1},\dots,o_{\underline{i}t_{\underline{i}}}$ be the
coefficient vectors corresponding to the missing and observed components of
$\hat{\bm{Z}}_{\underline{i}}$, respectively, and,
$\bm{c}^{m}_{\underline{i}}=\left(\bm{c}_{\underline{i}m_{\underline{i}1}}^{mT},\dots,\bm{c}_{\underline{i}m_{\underline{i}s_{\underline{i}}}}^{mT}\right)^{T}$
and
$\bm{c}_{\underline{i}}^{o}=\left(\bm{c}_{\underline{i}o_{\underline{i}1}}^{oT},\dots,\bm{c}_{\underline{i}o_{\underline{i}t_{\underline{i}}}}^{oT}\right)^{T}$.
Moreover, let indicate with $\hat{\bm{b}}_{lj}$, $l=1,\dots,L^{imp}$,
$j=1,\dots,p$, the coefficient vectors related to the basis expansions of the
components of the estimated principal components $\hat{\bm{\psi}}_{l}^{imp}$,
and with
$\hat{\bm{B}}=\left(\hat{\bm{b}}_{1},\dots,\hat{\bm{b}}_{L^{imp}}\right)$ the
matrix whose columns are
$\hat{\bm{b}}_{l}=\left(\hat{\bm{b}}_{l1}^{T},\dots,\hat{\bm{b}}_{lp}^{T}\right)^{T}$.
Then, the solution of the minimization problem in Equation (4) is

$\hat{\bm{c}}_{\underline{i}}^{m}=-\bm{C}_{mm}^{+}\bm{C}_{mo}\bm{c}_{\underline{i}}^{o}$
(5)

where $\bm{A}^{+}$ is the Moore-Penrose inverse of the matrix $\bm{A}$,
$\bm{C}_{mm}$ and $\bm{C}_{mo}$ are the matrices constructed by taking the
columns $m_{\underline{i}1},\dots,m_{\underline{i}s_{\underline{i}}}$ and
$o_{\underline{i}1},\dots,o_{\underline{i}t_{\underline{i}}}$ of the matrix
composed by the rows
$m_{\underline{i}1},\dots,m_{\underline{i}s_{\underline{i}}}$ of

$\bm{C}=\bm{W}\hat{\bm{B}}\hat{\bm{\Lambda}}^{-1}\hat{\bm{B}}^{T}\bm{W}$

with $\bm{W}$ the block-diagonal matrix defined in Section 2.1, and
$\hat{\bm{\Lambda}}$ the diagonal matrix whose diagonal entries are the
estimated eigenvalues $\hat{\lambda}_{l}^{imp}$, $l=1,\dots,L^{imp}$.
Moreover, to address the correlation bias issue typical of deterministic
imputation approaches (Little and Rubin, , 2019; Van Buuren, , 2018), we
propose to impute $\bm{c}_{\underline{i}}^{m}$ through
$\bm{c}_{\underline{i}}^{m,imp}=\left(\bm{c}_{\underline{i}1}^{m,impT},\dots,\bm{c}_{\underline{i}s_{\underline{i}}}^{m,impT}\right)^{T}$
as follows

$\bm{c}_{\underline{i}}^{m,imp}=\hat{\bm{c}}_{\underline{i}}^{m}+\bm{\varepsilon}_{i}$
(6)

where $\bm{\varepsilon}_{i}$ is a multivariate normal random variable with
mean zero and a residual covariance matrix robustly estimated from the
regression residuals of the coefficient vectors of the missing component on
those of the observed component for the observations in $S_{c}$ through the
model in Equation (5). Thus, the proposed RoFDI approach is a stochastic
imputation method (Little and Rubin, , 2019; Van Buuren, , 2018). Then, the
components of $\hat{\bm{Z}}_{\underline{i}}^{m}$ are imputed, for
$j=1,\dots,s_{\underline{i}}$, as

$\hat{Z}_{\underline{i}j}^{m,imp}(t)=\left(\bm{c}_{\underline{i}j}^{m,imp}\right)^{T}\bm{\phi}_{j}(t),\quad
t\in\mathcal{T},$

where $\bm{\phi}_{j}$ is the vector of basis function corresponding to the
$j$-th component of $\hat{\bm{Z}}$ (Section 2.1), and the imputed missing
components of $\bm{X}_{\underline{i}}$ are obtained by unstandardizing
$\hat{\bm{Z}}_{\underline{i}}^{m}$. Once the missing components of
$\bm{X}_{\underline{i}}$ are imputed, the whole observation is added to
$S_{c}$ and the next observation, which does not belong to $S_{c}$ and
corresponds to the smallest $s_{\underline{i}}$, is considered. Similarly to
Branden and Verboven,  (2009), if the cardinality of $S_{c}$ at the first
iteration is sufficiently large, we suggest to not update the RoMFPCA model
each time a new imputed observation is added to $S_{c}$ to avoid infeasible
time complexity of the RoFDI, otherwise, the RoMFPCA model could be updated
each time a given number of observations are added to $S_{c}$. Finally, to
take into account the increased noise due to single imputation, the proposed
RoFDI can be easily included in a multiple imputation framework (Van Buuren, ,
2018; Little and Rubin, , 2019), indeed, differently imputed datasets may be
obtained by performing several times the RoFDI due to the presence of the
stochastic component $\bm{\varepsilon}_{i}$ in Equation (6).

###  2.4 The Monitoring Strategy

The (IV) element of the proposed RoMFCC implementation relies on the
consolidated monitoring strategy for a multivariate functional quality
characteristic $\bm{X}$ based on the Hotelling’s $T^{2}$ and $SPE$ control
charts. The former assesses the stability of $\bm{X}$ in the finite
dimensional space spanned by the first principal components identified through
the RoMFPCA (Section 2.1), whereas, the latter monitors changes along
directions in the complement space. Specifically, the Hotelling’s $T^{2}$
statistic for $\bm{X}$ is defined as

$T^{2}=\sum_{l=1}^{L^{mon}}\frac{(\xi_{l}^{mon})^{2}}{\lambda_{l}^{mon}},$

where $\lambda_{l}^{mon}$ are the variances of the scores
$\xi_{l}^{mon}=\langle\bm{\psi}_{l}^{mon},\bm{Z}\rangle_{\mathbb{H}}$ where
$\bm{Z}$ is the vector of standardized variable of $\bm{X}$ and
$\bm{\psi}_{l}^{mon}$ are the corresponding principal components as defined in
Section 2.1. The number $L^{mon}$ is chosen such that the retained principal
components explain at least a given percentage $\delta^{mon}$ of the total
variability. The statistic $T^{2}$ is the standardized squared distance from
the centre of the orthogonal projection of $\bm{Z}$ onto the principal
component space spanned by
$\bm{\psi}_{1}^{mon},\dots,\bm{\psi}_{L^{mon}}^{mon}$. Whereas, the distance
between $\bm{Z}$ and its orthogonal projection onto the principal component
space is measured through the $SPE$ statistic, defined as

$SPE=||\bm{Z}-\hat{\bm{Z}}||_{\mathbb{H}}^{2},$

where $\hat{\bm{Z}}=\sum_{l=1}^{L^{mon}}\xi_{l}^{mon}\bm{\psi}_{l}^{mon}$.

Under the assumption of multivariate normality of $\xi_{l}^{mon}$, which is
approximately true by the central limit theorem (Nomikos and MacGregor, ,
1995), the control limits of the Hotelling’s $T^{2}$ control chart can be
obtained by considering the $(1-\alpha^{*})$ quantiles of a chi-squared
distribution with $L^{mon}$ degrees of freedom (Johnson et al., , 2014).
Whereas, the control limits for $SPE$ control chart can be computed by using
the following equation (Jackson and Mudholkar, , 1979)

$CL_{SPE,\alpha^{*}}=\theta_{1}\left[\frac{c_{\alpha^{*}}\sqrt{2\theta_{2}h_{0}^{2}}}{\theta_{1}}+1+\frac{\theta_{2}h_{0}(h_{0}-1)}{\theta_{1}^{2}}\right]^{1/h_{0}}$

where $c_{\alpha^{*}}$ is the normal deviate corresponding to the upper
$(1-\alpha^{*})$ quantile, $h_{0}=1-2\theta_{1}\theta_{3}/3\theta_{2}^{2}$,
$\theta_{j}=\sum_{l=L^{mon}+1}^{\infty}(\lambda_{l}^{mon})^{j}$, $j=1,2,3$.
Note that, to control the family wise error rate (FWER), $\alpha^{*}$ should
be chosen appropriately. We propose to use the Šidák correction
$\alpha^{*}=1-\left(1-\alpha\right)^{1/2}$ (Šidák, , 1967), where $\alpha$ is
the overall type I error probability.

###  2.5 The Proposed Method

The proposed RoMFCC implementation collects all the elements introduced in the
previous sections for the Phase II monitoring strategy where a set of Phase I
observations, which can be contaminated with both functional casewise and
cellwise outliers, is used for the design of the control chart. Both phases
are outlined in the scheme of Figure 2 and detailed in the following sections.

Figure 2: Scheme of the RoMFCC approach.

####  2.5.1 Phase I

Let $\bm{X}_{i}$, $i=1,\dots,n$, the Phase I random sample of the multivariate
functional quality characteristic $\bm{X}$, used to characterize the normal
operating conditions of the process. It can be contaminated with both
functional casewise and cellwise outliers.

In the filtering step, functional cellwise outliers are identified through the
FUF described in Section 2.2, and, then, replaced by missing components. Note
that, if cellwise outliers are identified for each component of a given
observation, then, that observation is removed from the sample because its
imputation does not provide any additional information for the analysis. In
the imputation step, missing components are imputed through the RoFDI method
presented in Section 2.3. Once the imputed Phase I sample is obtained, it used
to to estimate the RoMFPCA model and perform the dimensionality reduction step
as described in Section 2.1. The Hotelling’s $T^{2}$ and $SPE$ statistics are,
then, computed for each observation in the Phase I sample after the imputation
step. Specifically, the values $T^{2}_{i}$ and $SPE_{i}$ of the statistics are
computed as described in Section 2.4 by considering the estimated RoMFPCA
model obtained in the dimensionality reduction step. Finally, control limits
for the Hotelling’s $T^{2}$ and $SPE$ control charts are obtained as described
in Section 2.4. Note that the parameters $\theta_{j}$ to estimate
$CL_{SPE,\alpha^{*}}$ are approximated by considering a finite summation to
the maximum number of estimable principal components, which is finite for a
sample of $n$ observations. If a multiple imputation strategy is employed by
performing the imputation step several times, then the multiple estimated
RoMFPCA models could be combined by averaging the robustly estimated
covariance functions as suggested in Van Ginkel et al.,  (2007).

Note that, when the sample size $n$ is small compared to the number of process
variables, undesirable effect upon the performance of the RoMFCC could arise
(Ramaker et al., , 2004; Kruger and Xie, , 2012). To reduce possible
overfitting issues and, thus, increase the monitoring performance of the
RoMFCC, a reference sample of Phase I observations, referred to as tuning set,
could be considered, which is different from the one used in the previous
steps, referred to as training set. Specifically, on the line of Kruger and
Xie,  (2012), Chapter 6.4, the tuning set is passed through the filtering and
imputation steps and, then, it is projected on the RoMFPCA model, estimated
through the training set observations, to robustly estimate the distribution
of the resulting scores. Finally, Hotelling’s $T^{2}$ and $SPE$ statistics and
control charts limits are calculated by taking into account the estimated
distribution of the tuning set scores.

####  2.5.2 Phase II

In the actual monitoring phase (Phase II), a new observation $\bm{X}_{new}$ is
projected on the RoMFPCA model to compute the values of $T^{2}_{new}$ and
$SPE_{new}$ statistics accordingly to the score distribution identified in
Phase I. An alarm signal is issued if at least one between $T^{2}_{new}$ and
$SPE_{new}$ violates the control limits.

##  3 Simulation Study

The overall performance of the proposed RoMFCC is evaluated by means of an
extensive Monte Carlo simulation study. The aim of this simulation is to
assess the performance of the RoMFCC in identifying mean shifts of the
multivariate functional quality characteristic when the Phase I sample is
contaminated with both functional cellwise and casewise outliers. The data
generation process (detailed in Supplementary Material A) is inspired by the
real-case study in Section 4 and mimics typical behaviours of DRCs in a RSW
process. Specifically, it considers a multivariate functional quality
characteristic with $p=10$ components. Two main scenarios are considered that
are characterized by different Phase I sample contamination. Specifically, the
Phase I sample is contaminated by functional cellwise outliers in Scenario 1
and by functional casewise outliers in Scenario 2, with a contamination
probability equal to 0.05 in both cases. In the Supplementary Material B,
additional results obtained by considering contamination probability equal to
0.1 are shown. For each scenario, two contamination models, referred to as Out
E and Out P, with three increasing contamination levels, referred to as C1, C2
and C3, are considered. The former mimics a splash weld (expulsion), caused by
excessive welding current, while the latter resembles phase shift of the peak
time caused by an increased electrode force (Xia et al., , 2019). Moreover, we
consider also a scenario, referred to as Scenario 0, representing settings
where the Phase I sample is not contaminated. To generate the Phase II sample,
two types of OC conditions, referred to as OC E and OC P, are considered that
are generated analogously to the two contamination models Out E and Out P,
respectively, at 4 different severity levels $SL=\\{1,2,3,4\\}$.

The proposed RoMFCC implementation is compared with several natural competing
approaches. The first approaches are control charts for multivariate scalar
data, built on the average value of each component of the multivariate
functional data. Among them, we consider the multivariate classical
Hotelling’s $T^{2}$ control chart, referred to as M, the multivariate
iterative variant, referred to as Miter, where outliers detected by the
control chart in Phase I are iteratively removed until all data are assumed to
be IC, and the multivariate robust control chart proposed by Chenouri et al.,
(2009), referred to as MRo. We further consider also two approaches recently
appeared in the profile monitoring literature, i.e., the multivariate
functional control charts, referred to as MFCC, proposed by Capezza et al.,
(2020, 2022), and the multivariate iterative functional control charts
variant, referred to as MFCCiter, where outliers detected by the control chart
in Phase I are iteratively removed until all data are assumed to be IC. The
RoMFCC is implemented as described in Section 2 with
$\delta_{fil}=\delta_{imp}=0.999$ and $\delta_{mon}=0.7$, and to take into
account the increased noise due to single imputation, 5 differently imputed
datasets are generated through RoFDI. While data are observed through noisy
discrete values, each component of the generated quality characteristic
observations is obtained by considering the approximation in Equation (1) with
$K=10$ cubic B-splines estimated through the spline smoothing approach based
on a roughness penalty on the second derivative (Ramsay and Silverman, ,
2005). For each scenario, contamination model, contamination level, OC
condition and severity level, 50 simulation runs are performed. Each run
considers a Phase I sample of 4000 observations, where for MFCC, MFCCIter and
RoMFCC, 1000 are used as training set, and the remaining 3000 are used as
tuning set. The Phase II sample is composed of 4000 i.i.d. observations. The
RoMFCC and the competing methods performances are assessed by means of true
detection rate (TDR), which is the proportion of points outside the control
limits whilst the process is OC, and the false alarm rate (FAR), which is the
proportion of points outside the control limits whilst the process is IC. The
FAR should be as similar as possible to the overall type I error probability
$\alpha$ considered to obtain the control limits and set equal to 0.05,
whereas the TDR should be as close to one as possible.

Figure 3-5 display for Scenario 0, Scenario 1 and Scenario 2, respectively, as
a function of the severity level $SL$, the mean FAR ($SL=0$) or TDR ($SL\neq
0$) for each OC condition OC E and OC P, contamination level C1, C2 and C3 and
contamination model Out E and Out P.

Figure 3: Mean FAR ($SL=0$) or TDR ($SL\neq 0$) achieved by M, Miter, MRo, MFCC, MFCCiter and RoMFCC for each OC condition (OC E and OC P) as a function of the severity level $SL$ in Scenario 0. OC E |  OC P
---|---
|
Figure 4: Mean FAR ($SL=0$) or TDR ($SL\neq 0$) achieved by M, Miter, MRo, MFCC, MFCCiter and RoMFCC for each contamination level (C1, C2 and C3), OC condition (OC E and OC P) as a function of the severity level $SL$ with contamination model Out E and Out P in Scenario 1. |  Out E |  Out P
---|---|---
|  OC E |  OC P |  OC E |  OC P
C1 |  |  |  |
C2 |  |  |  |
C3 |  |  |  |
Figure 5: Mean FAR ($SL=0$) or TDR ($SL\neq 0$) achieved by M, Miter, MRo, MFCC, MFCCiter and RoMFCC for each contamination level (C1, C2 and C3), OC condition (OC E and OC P) as a function of the severity level $SL$ with contamination model Out E and Out P in Scenario 2. |  Out E |  Out P
---|---|---
|  OC E |  OC P |  OC E |  OC P
C1 |  |  |  |
C2 |  |  |  |
C3 |  |  |  |

When the Phase I sample is not contaminated by outliers, Figure 3 shows that
all the approaches that take into account the functional nature of the data,
i.e., MFCC, MFCCiter, RoMFCC, achieve the same performance for both OC
conditions. Although this setting should be not favourable to approaches
specifically designed to deal with outliers, MFCCiter and RoMFCC perform equal
to MFCC. The non-functional approaches, i.e., M, Miter, MRo, show worse
performance than the functional counterparts, and there is no significant
performance difference among them as well. Figure 4 show the results for
Scenario 1, where the Phase I sample is contaminated by cellwise outliers. The
proposed RoMFCC largely outperforms the competing methods for each
contmination model, contamination level and OC condition. As expected, as the
contamination level increases the differences in performance between the
RoMFCC and the competing methods increase as well. Indeed, the performance of
RoMFCC is almost insensible to the contamination levels as well as the
contamination models, differently from the competing methods whose performance
decreases as the contamination level increases and the contamination model
changes. Moreover, the MFCCiter, which is representative of the baseline
method in this setting, only slightly improves the performance of the MFCC.
This is probably due to the masking effect that prevents the MFCCiter to
iteratively identify functional cellwise outlier in the Phase I sample and,
thus, makes it equivalent to MFCC. The performance of the non-functional
methods are totally unsatisfactory because they are not able to both capture
the functional nature of the data and successfully deal with functional
cellwise outliers. Specifically, M is the worst method overall readily
followed by Miter and MRo. Furthermore, as the contamination level increases,
the competing functional methods lose their favourable performance and tend to
achieve performance comparable to the non-functional approaches. This shows
that if outliers are not properly dealt with, terrible performance could arise
independently of the suitability of the methods considered.

From Figure 5, also in Scenario 2, RoMFCC is clearly the best method. However,
in this scenario the difference in performance between the RoMFCC and the
competing methods is less pronounced than in Scenario 1. This is expected
because the Phase I sample is now contaminated by functional casewise
outliers, which is the only contamination type against which the competing
methods are able to provide a certain degree of robustness. Note that, the
performance of RoMFCC are almost unaffected by contamination in the Phase I
sample also in this scenario, which proves the ability of the proposed method
to deal with both functional cellwise and casewise outliers.

##  4 Real-Case Study

To demonstrate the potential of the proposed RoMFCC in practical situations, a
real-case study in the automotive industry is presented henceforth. As
introduced in Section 1, it addresses the issue of monitoring the quality of
the RSW process, which guarantees the structural integrity and solidity of
welded assemblies in each vehicle (Martín et al., , 2014). The RSW process
(Zhang and Senkara, , 2011) is an autogenous welding process in which two
overlapping conventional steel galvanized sheets are joint together, without
the use of any filler material. Joints are formed by applying pressure to the
weld area from two opposite sides by means of two copper electrodes. Voltage
applied to the electrodes generates a current flowing between them through the
material. The electrical current flows because the resistance offered by
metals causes significant heat generation (Joule effect) that increases the
metal temperature at the faying surfaces of the work pieces up to the melting
point. Finally, due to the mechanical pressure of the electrodes, the molten
metal of the jointed metal sheets cools and solidifies, forming the so-called
weld nugget (Raoelison et al., , 2012). To monitor the RSW process, the modern
automotive Industry 4.0 framework allows the automatic acquisition of a large
volume of process parameters. The DRC is considered the most important of
these parameters to describe the quality of the RSW process. Further details
on how the typical behaviour of a DRC is related to the physical and
metallurgical development of a spot weld are provided by Capezza et al., 2021b
.

Data analyzed in this study are courtesy of Centro Ricerche Fiat and are
recorded at the Mirafiori Factory during lab tests and are acquired during RSW
processes made on the body of the Fiat 500BEV. A body is characterized by a
large number of spot welds with different characteristics, e.g, the thickness
and material of the sheets to be joined together and the welding time. In this
paper, we focus on monitoring a set of ten spot welds made on the body by one
of the welding machines. Therefore, for each sub-assembly the multivariate
functional quality characteristic is a vector of the ten DRCs, corresponding
to the the second welding pulse of ten spot welds normalized on the time
domain $[0,1]$, for a total number of assemblies equal to 1839. Moreover,
resistance measurements were collected at a regular grid of points equally
spaced by 1 ms.

The RSW process quality is directly affected by electrode wear since the
increase in weld numbers leads to changed electrical, thermal and mechanical
contact conditions at electrode and sheet interfaces (Manladan et al., ,
2017). Thus, to take into account the wear issue, electrodes go through
periodical renovations. In this setting, a paramount issue refers to the swift
identification of DRCs mean shifts caused by electrode wear, which could be
considered as a criterion for electrode life termination and guide the
electrode renovation strategy.

In the light of this, the 919 multivariate profiles corresponding to spot
welds made immediately before electrode renewal are used to form the Phase I
sample, whereas, the remaining 920 observations are used in Phase II to
evaluate the proposed chart performance. We expect that the mean shift of the
Phase II DRCs caused by electrode wear should be effectively captured by the
proposed control chart. The RoMFCC is implemented as in the simulation study
in Section 3 with the training and tuning sets, each composed by 460 Phase I
observations, randomly selected without remittance. As shown in Figure 1 of
Section 1, data are plausibly contaminated by several outliers. This is
further confirmed by Figure 6, which shows the boxplot of the functional
distance square roots $\sqrt{D_{i,fil}}$ (Equation (3)) obtained from the FUF
applied on the training set. Some components clearly show the presence of
functional cellwise outliers that are possibly arranged in groups, while other
components seem less severely contaminated.

Figure 6: Boxplot of the functional distance square roots $\sqrt{D_{i,fil}}$
(Equation (3)) obtained from the FUF applied on the training set.

Figure 7 shows the application of the proposed RoMFCC.

Figure 7: Hotelling’s $T^{2}$ and $SPE$ control charts for the RoMFCC in the
real-case study. The vertical line separates the monitoring statistics
calculated for the tuning set, on the left, and the Phase II data set on the
right, while the horizontal lines define the control limits.

The vertical line separates the monitoring statistics calculated for the
tuning set, on the left, and the Phase II data set on the right, while the
horizontal lines define the control limits. Note that, a significant number of
tuning set observations are signaled as OC, highlighted in red in Figure 7.
This is expected because these points may include functional casewise outlier
not filtered out by the FUF. In the monitoring phase, many points are signaled
as OCs by the RoMFCC. In particular, the RoMFCC signals 72.3% of the
observations in the Phase II data set as OC. This shows that the proposed
method is particularly sensible to mean shifts caused by an increased
electrode wear.

Finally, the proposed method is compared with the competing methods presented
in the simulation study in Section 3. Table 1 shows the estimated TDR values
$\widehat{TDR}$ on the Phase II sample for all the considered competing
methods. Similarly to Centofanti et al., 2021b , the uncertainty of
$\widehat{TDR}$ is quantified through a bootstrap analysis (Efron and
Tibshirani, , 1994). Table 1 reports the mean $\overline{TDR}$ of the
empirical bootstrap distribution of $\widehat{TDR}$, and the corresponding
bootstrap 95% confidence interval (CI) for each monitoring method.

| $\widehat{TDR}$ | $\overline{TDR}$ | CI
---|---|---|---
M | 0.336 | 0.335 | [0.305,0.368]
Miter | 0.462 | 0.461 | [0.428,0.496]
MRo | 0.513 | 0.512 | [0.481,0.547]
MFCC | 0.541 | 0.541 | [0.511,0.574]
MFCCiter | 0.632 | 0.632 | [0.595,0.664]
RoMFCC | 0.723 | 0.723 | [0.695,0.753]

Table 1: Estimated TDR values $\hat{TDR}$ on the Phase II sample, mean
$\overline{T}DR$ of the empirical bootstrap distribution of $\hat{TDR}$, and
the corresponding bootstrap 95% confidence interval (CI) for each monitoring
method in the real-case study.

The bootstrap analysis shows that the RoMFCC outperforms the competing control
charts, indeed bootstrap 95% confidence intervals are strictly above those of
all considered monitoring approaches. As in Section 3, non-functional
approaches, i.e., M, Miter, MRo, show worse performance than the functional
counterparts because they are not able to satisfactorily capture the
functional nature of the data and robust approaches always improve the non-
robust ones. Therefore, the proposed RoMFCC stands out as the best method to
promptly identify OC conditions in the RWS process caused by an increased
electrode wear with a Phase I sample contaminated by functional outliers.

##  5 Conclusions

In this paper, we propose a new robust framework for the statistical process
monitoring of multivariate functional data, referred to as robust multivariate
functional control charts (RoMFCC). The RoMFCC is designed to assess the
presence of assignable causes of variation while being robust to both
functional casewise and cellwise outliers. The proposed method is suitable for
those industrial processes where many functional variables are available and
occasional outliers are produced, such as anomalies in the data acquisition
and data collected during a fault in the process. Specifically, the RoMFCC
framework is based on four main elements, i.e. a functional univariate filter
to identify functional cellwise outliers, to be replaced by missing values, a
robust functional data imputation of these missing values, a casewise robust
dimensionality reduction based on ROBPCA and a monitoring strategy based on
the Hotelling’s $T^{2}$ and $SPE$ control charts. These elements are combined
in a Phase II monitoring strategy where a set of Phase I observations, which
can be contaminated with both functional casewise and cellwise outliers, is
used for the design of the control chart.

To the best of the authors’ knowledge, the RoMFCC framework is the first
monitoring scheme that is able to monitor a multivariate functional quality
characteristic while being robust to functional casewise and cellwise
outliers. Indeed, methods already present in the literature either apply
robust approaches to multivariate scalar features extracted from the profiles
or use diagnostic approaches on the multivariate functional data to
iteratively remove outliers. However, the former are not able to capture the
functional nature of the data, while the latter are not able to deal with
functional cellwise outliers.

The performance of the RoMFCC framework is assessed through an extensive Monte
Carlo simulation study where it is compared with several competing monitoring
methods for multivariate scalar data and multivariate functional data. The
ability of the proposed method to estimate the distribution of the data
without removing observations while being robust to both functional casewise
and cellwise outliers allows the RoMFCC to outperform the competitors in all
the considered scenarios. Lastly, the practical applicability of the proposed
method is illustrated through a motivating real-case study, which addresses
the issue of monitoring the quality of a resistance spot-welding process in
the automotive industry. Also in this case, the RoMFCC shows better
performance than the competitors in the identification of out-of-control
condition of the dynamic resistance curves.

## Supplementary Materials

The Supplementary Materials contain additional details about the data
generation process in the simulation study (A), additional simulation results
(B), as well as the R code to reproduce graphics and results over competing
methods in the simulation study.

## Acknowledgments

The present work was developed within the activities of the project
ARS01_00861 “Integrated collaborative systems for smart factory - ICOSAF”
coordinated by CRF (Centro Ricerche Fiat Scpa - www.crf.it) and financially
supported by MIUR (Ministero dell’Istruzione, dell’Università e della
Ricerca).

## References

  * Aerts and Wilms, (2017) Aerts, S. and Wilms, I. (2017).  Cellwise robust regularized discriminant analysis.  Statistical Analysis and Data Mining: The ASA Data Science Journal, 10(6):436–447.
  * Agostinelli et al., (2015) Agostinelli, C., Leung, A., Yohai, V. J., and Zamar, R. H. (2015).  Robust estimation of multivariate location and scatter in the presence of cellwise and casewise contamination.  Test, 24(3):441–461.
  * Alemán-Gómez et al., (2022) Alemán-Gómez, Y., Arribas-Gil, A., Desco, M., Elías, A., and Romo, J. (2022).  Depthgram: Visualizing outliers in high-dimensional functional data with application to fmri data exploration.  Statistics in Medicine.
  * Alfaro and Ortega, (2009) Alfaro, J. and Ortega, J. F. (2009).  A comparison of robust alternatives to hotelling’s $t^{2}$ control chart.  Journal of Applied Statistics, 36(12):1385–1396.
  * Alqallaf et al., (2009) Alqallaf, F., Van Aelst, S., Yohai, V. J., and Zamar, R. H. (2009).  Propagation of outliers in multivariate data.  The Annals of Statistics, pages 311–331.
  * Arribas-Gil and Romo, (2014) Arribas-Gil, A. and Romo, J. (2014).  Shape outlier detection and visualization for functional data: the outliergram.  Biostatistics, 15(4):603–619.
  * Boente and Salibián-Barrera, (2021) Boente, G. and Salibián-Barrera, M. (2021).  Robust functional principal components for sparse longitudinal data.  Metron, 79(2):159–188.
  * Branden and Verboven, (2009) Branden, K. V. and Verboven, S. (2009).  Robust data imputation.  Computational Biology and Chemistry, 33(1):7–13.
  * Cabana and Lillo, (2021) Cabana, E. and Lillo, R. E. (2021).  Robust multivariate control chart based on shrinkage for individual observations.  Journal of Quality Technology, pages 1–26.
  * (10) Capezza, C., Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2021a).  Functional regression control chart for monitoring ship co2 emissions.  Quality and Reliability Engineering International, 38(3):1519–1537.
  * Capezza et al., (2022) Capezza, C., Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2022).  funcharts: Control charts for multivariate functional data in r.
  * (12) Capezza, C., Centofanti, F., Lepore, A., and Palumbo, B. (2021b).  Functional clustering methods for resistance spot welding process data in the automotive industry.  Applied Stochastic Models in Business and Industry, 37(5):908–925.
  * Capezza et al., (2020) Capezza, C., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2020).  Control charts for monitoring ship operating conditions and co2 emissions based on scalar-on-function regression.  Applied Stochastic Models in Business and Industry, 36(3):477–500.
  * (14) Centofanti, F., Colosimo, B. M., Grasso, M. L., Menafoglio, A., Palumbo, B., and Vantini, S. (2021a).  Robust functional anova with application to additive manufacturing.
  * (15) Centofanti, F., Lepore, A., Menafoglio, A., Palumbo, B., and Vantini, S. (2021b).  Functional regression control chart.  Technometrics, 63(3):281–294.
  * Chenouri et al., (2009) Chenouri, S., Steiner, S. H., and Variyath, A. M. (2009).  A multivariate robust control chart for individual observations.  Journal of Quality Technology, 41(3):259–271.
  * Chiou et al., (2014) Chiou, J.-M., Chen, Y.-T., and Yang, Y.-F. (2014).  Multivariate functional principal component analysis: A normalization approach.  Statistica Sinica, pages 1571–1596.
  * Cuesta-Albertos and Fraiman, (2006) Cuesta-Albertos, J. A. and Fraiman, R. (2006).  Impartial trimmed means for functional data.  DIMACS Series in Discrete Mathematics and Theoretical Computer Science, 72:121.
  * Cuevas and Fraiman, (2009) Cuevas, A. and Fraiman, R. (2009).  On depth measures and dual statistics. a methodology for dealing with general data.  Journal of Multivariate Analysis, 100(4):753–766.
  * Dai and Genton, (2018) Dai, W. and Genton, M. G. (2018).  Multivariate functional data visualization and outlier detection.  Journal of Computational and Graphical Statistics, 27(4):923–934.
  * Dickinson et al., (1980) Dickinson, D., Franklin, J., Stanya, A., et al. (1980).  Characterization of spot welding behavior by dynamic electrical parameter monitoring.  Welding Journal, 59(6):170.
  * Efron and Tibshirani, (1994) Efron, B. and Tibshirani, R. J. (1994).  An introduction to the bootstrap.  CRC press.
  * Febrero et al., (2008) Febrero, M., Galeano, P., and González-Manteiga, W. (2008).  Outlier detection in functional data by depth measures, with application to identify abnormal nox levels.  Environmetrics: The official journal of the International Environmetrics Society, 19(4):331–345.
  * Filzmoser et al., (2020) Filzmoser, P., Höppner, S., Ortner, I., Serneels, S., and Verdonck, T. (2020).  Cellwise robust m regression.  Computational Statistics & Data Analysis, 147:106944.
  * Fraiman and Muniz, (2001) Fraiman, R. and Muniz, G. (2001).  Trimmed means for functional data.  Test, 10(2):419–440.
  * Grasso et al., (2016) Grasso, M., Menafoglio, A., Colosimo, B. M., and Secchi, P. (2016).  Using curve-registration information for profile monitoring.  Journal of Quality Technology, 48(2):99–127.
  * Happ and Greven, (2018) Happ, C. and Greven, S. (2018).  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  Journal of the American Statistical Association, 113(522):649–659.
  * Hubert et al., (2015) Hubert, M., Rousseeuw, P. J., and Segaert, P. (2015).  Multivariate functional outlier detection.  Statistical Methods & Applications, 24(2):177–202.
  * Hubert et al., (2005) Hubert, M., Rousseeuw, P. J., and Vanden Branden, K. (2005).  Robpca: a new approach to robust principal component analysis.  Technometrics, 47(1):64–79.
  * Hyndman and Shang, (2010) Hyndman, R. J. and Shang, H. L. (2010).  Rainbow plots, bagplots, and boxplots for functional data.  Journal of Computational and Graphical Statistics, 19(1):29–45.
  * Hyndman and Ullah, (2007) Hyndman, R. J. and Ullah, M. S. (2007).  Robust forecasting of mortality and fertility rates: A functional data approach.  Computational Statistics & Data Analysis, 51(10):4942–4956.
  * Jackson and Mudholkar, (1979) Jackson, J. E. and Mudholkar, G. S. (1979).  Control procedures for residuals associated with principal component analysis.  Technometrics, 21(3):341–349.
  * Jensen et al., (2007) Jensen, W. A., Birch, J. B., and Woodall, W. H. (2007).  High breakdown estimation methods for phase i multivariate control charts.  Quality and Reliability Engineering International, 23(5):615–629.
  * Johnson et al., (2014) Johnson, R. A., Wichern, D. W., et al. (2014).  Applied multivariate statistical analysis, volume 6.  Pearson London, UK:.
  * Jolliffe, (2011) Jolliffe, I. (2011).  Principal component analysis.  Springer.
  * Kokoszka and Reimherr, (2017) Kokoszka, P. and Reimherr, M. (2017).  Introduction to functional data analysis.  Chapman and Hall/CRC.
  * Kordestani et al., (2020) Kordestani, M., Hassanvand, F., Samimi, Y., and Shahriari, H. (2020).  Monitoring multivariate simple linear profiles using robust estimators.  Communications in Statistics-Theory and Methods, 49(12):2964–2989.
  * Kruger and Xie, (2012) Kruger, U. and Xie, L. (2012).  Statistical monitoring of complex multivatiate processes: with applications in industrial process control.  John Wiley & Sons.
  * Lee et al., (2013) Lee, S., Shin, H., and Billor, N. (2013).  M-type smoothing spline estimators for principal functions.  Computational Statistics & Data Analysis, 66:89–100.
  * Leung et al., (2016) Leung, A., Zhang, H., and Zamar, R. (2016).  Robust regression estimation and inference in the presence of cellwise and casewise contamination.  Computational Statistics & Data Analysis, 99:1–11.
  * Little and Rubin, (2019) Little, R. J. and Rubin, D. B. (2019).  Statistical analysis with missing data, volume 793.  John Wiley & Sons.
  * Locantore et al., (1999) Locantore, N., Marron, J., Simpson, D., Tripoli, N., Zhang, J., Cohen, K., Boente, G., Fraiman, R., Brumback, B., Croux, C., et al. (1999).  Robust principal component analysis for functional data.  Test, 8(1):1–73.
  * López-Pintado and Romo, (2011) López-Pintado, S. and Romo, J. (2011).  A half-region depth for functional data.  Computational Statistics & Data Analysis, 55(4):1679–1695.
  * Manladan et al., (2017) Manladan, S., Yusof, F., Ramesh, S., Fadzil, M., Luo, Z., and Ao, S. (2017).  A review on resistance spot welding of aluminum alloys.  The International Journal of Advanced Manufacturing Technology, 90(1):605–634.
  * Maronna et al., (2019) Maronna, R. A., Martin, R. D., Yohai, V. J., and Salibián-Barrera, M. (2019).  Robust statistics: theory and methods (with R).  John Wiley & Sons.
  * Martín et al., (2014) Martín, Ó., Pereda, M., Santos, J. I., and Galán, J. M. (2014).  Assessment of resistance spot welding quality based on ultrasonic testing and tree-based techniques.  Journal of Materials Processing Technology, 214(11):2478–2487.
  * Menafoglio et al., (2018) Menafoglio, A., Grasso, M., Secchi, P., and Colosimo, B. M. (2018).  Profile monitoring of probability density functions via simplicial functional pca with application to image data.  Technometrics, 60(4):497–510.
  * Moheghi et al., (2020) Moheghi, H., Noorossana, R., and Ahmadi, O. (2020).  Phase i and phase ii analysis of linear profile monitoring using robust estimators.  Communications in Statistics-Theory and Methods, pages 1–18.
  * Montgomery, (2012) Montgomery, D. C. (2012).  Introduction to Statistical Quality Control.  Wiley.
  * Nomikos and MacGregor, (1995) Nomikos, P. and MacGregor, J. F. (1995).  Multivariate spc charts for monitoring batch processes.  Technometrics, 37(1):41–59.
  * Noorossana et al., (2011) Noorossana, R., Saghaei, A., and Amiri, A. (2011).  Statistical analysis of profile monitoring, volume 865.  John Wiley & Sons.
  * R Core Team, (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Ramaker et al., (2004) Ramaker, H.-J., van Sprang, E. N., Westerhuis, J. A., and Smilde, A. K. (2004).  The effect of the size of the training set and number of principal components on the false alarm rate in statistical process monitoring.  Chemometrics and intelligent laboratory systems, 73(2):181–187.
  * Ramsay and Silverman, (2005) Ramsay, J. O. and Silverman, B. W. (2005).  Functional Data Analysis.  Springer.
  * Raoelison et al., (2012) Raoelison, R., Fuentes, A., Rogeon, P., Carre, P., Loulou, T., Carron, D., and Dechalotte, F. (2012).  Contact conditions on nugget development during resistance spot welding of zn coated steel sheets using rounded tip electrodes.  Journal of Materials Processing Technology, 212(8):1663–1669.
  * Rousseeuw, (1984) Rousseeuw, P. J. (1984).  Least median of squares regression.  Journal of the American statistical association, 79(388):871–880.
  * Rousseeuw and Bossche, (2018) Rousseeuw, P. J. and Bossche, W. V. D. (2018).  Detecting deviating data cells.  Technometrics, 60(2):135–145.
  * Sawant et al., (2012) Sawant, P., Billor, N., and Shin, H. (2012).  Functional outlier detection with robust functional principal component analysis.  Computational Statistics, 27(1):83–102.
  * Šidák, (1967) Šidák, Z. (1967).  Rectangular confidence regions for the means of multivariate normal distributions.  Journal of the American Statistical Association, 62(318):626–633.
  * Sinova et al., (2018) Sinova, B., Gonzalez-Rodriguez, G., and Van Aelst, S. (2018).  M-estimators of location for functional data.  Bernoulli, 24(3):2328–2357.
  * Tarr et al., (2016) Tarr, G., Müller, S., and Weber, N. C. (2016).  Robust estimation of precision matrices under cellwise contamination.  Computational Statistics & Data Analysis, 93:404–420.
  * Van Buuren, (2018) Van Buuren, S. (2018).  Flexible imputation of missing data.  CRC press.
  * Van Ginkel et al., (2007) Van Ginkel, J. R., Van der Ark, L. A., Sijtsma, K., and Vermunt, J. K. (2007).  Two-way imputation: A bayesian method for estimating missing scores in tests and questionnaires, and an accurate approximation.  Computational Statistics & Data Analysis, 51(8):4013–4027.
  * Vargas, (2003) Vargas, N. J. A. (2003).  Robust estimation in multivariate control charts for individual observations.  Journal of Quality Technology, 35(4):367–376.
  * Xia et al., (2019) Xia, Y.-J., Su, Z.-W., Li, Y.-B., Zhou, L., and Shen, Y. (2019).  Online quantitative evaluation of expulsion in resistance spot welding.  Journal of Manufacturing Processes, 46:34–43.
  * Zhang and Senkara, (2011) Zhang, H. and Senkara, J. (2011).  Resistance welding: fundamentals and applications.  CRC press.
  * Zhou and Cai, (2014) Zhou, K. and Cai, L. (2014).  Study on effect of electrode force on resistance spot welding process.  Journal of applied physics, 116(8):084902.