Yuanyong Liang

Arnold Galindo

Chike Yuan

High Performance Interconnect Technologies for Supercomputing

Keywords:

High-performance computing (HPC) is contributing significantly to scientific progression thanks to its abundant computing power. With recent breakthroughs in natural language processing (NLP), HPC has become indispensable in driving their success, owing to its capacity to train large language models (LLMs). HPC systems not only demonstrate their importance in scientific computing but also in machine learning and generative artificial intelligence. This unprecedented computing power is firstly due to the development of newer computing resources: CPUs, GPUs, FPGAs, etc. However, the progress in interconnect networking, which links hundreds of thousands of compute nodes, should not be overlooked. This survey investigates current popular interconnect topologies driving the most powerful supercomputers.

INTRODUCTION

High-Performance Computing (HPC) interconnect topology refers to the arrangement and organization of communication paths in a high-performance computing system. The topology plays a crucial role in determining how different nodes (computers, processors, or other computing elements) within an HPC system communicate with each other. This communication is essential for the efficient processing and transfer of data, which is at the heart of HPC operations [START_REF] Abts | High performance datacenter networks: Architectures, algorithms, and opportunities[END_REF].

The effectiveness of an HPC interconnect topology is measured by several factors, including its bandwidth (the rate at which data is transferred), latency (the time delay in transferring data), scalability (how well the system performs as more nodes are added), and fault tolerance (the system's ability to continue operating in the event of a failure) [START_REF] James | Principles and practices of interconnection networks[END_REF].

There are several types of HPC interconnect topologies, such as fat-tree [START_REF] Charles E Leiserson | Fat-trees: Universal networks for hardware-efficient supercomputing[END_REF], dragonfly [START_REF] Kim | Technology-driven, highly-scalable dragonfly topology[END_REF][START_REF] Kang | Study of workload interference with intelligent routing on dragonfly[END_REF], and HyperX [START_REF] Ho Ahn | Hyperx: topology, routing, and packaging of efficient large-scale networks[END_REF]. Each topology is designed to serve special purposes with its own set of advantages and disadvantages. The choice of topology depends on the specific requirements and goals of the HPC system, such as the type of computations it will perform, the size of the datasets, budget constraints, and the desired balance between speed, cost, and reliability.

In recent years, there has been significant innovation in HPC interconnect technologies, driven by the growing demands for faster and more efficient computing capabilities, especially in fields like climate modeling, genomic analysis, and artificial intelligence [START_REF] Achiam | [END_REF][START_REF] Wang | Union: An automatic workload manager for accelerating network simulation[END_REF][START_REF] Van | Gromacs: fast, flexible, and free[END_REF][START_REF] James C Phillips | Scalable molecular dynamics with namd[END_REF]. These advances continue to push the boundaries of what's possible in high-performance computing. * All authors contributed equally to this work.

FAT-TREE TOPOLOGY

Fat-tree topology is a hierarchical network design that mimics a tree structure, where each layer in the hierarchy is connected upward to the next layer until it reaches the root [START_REF] Charles E Leiserson | Fat-trees: Universal networks for hardware-efficient supercomputing[END_REF][START_REF] Petrini | k-ary n-trees: High performance networks for massively parallel architectures[END_REF]. This topology is widely used in large-scale and high-performance computing environments, such as data centers and supercomputers, due to its scalability, high bandwidth, and fault tolerance capabilities. In the following sections, we will discuss the architecture of fat-tree topology, its advantages, challenges, and applications.

Architecture

The architecture of a fat-tree topology is designed to overcome the limitations of traditional tree topologies, especially in terms of bandwidth and scalability. In a fat-tree, the concept of a "fat" tree is realized by increasing the number of links and the bandwidth as one moves up the hierarchy. This is in contrast to a conventional tree where all links typically have the same bandwidth. As depicted in Figure 1, the bottom layer of a fat-tree consists of edge switches, which connect to nodes such as servers or computers. The next layer consists of aggregation switches, which connect multiple edge switches. Finally, at the top of the hierarchy, core switches interconnect the aggregation switches. The key feature of this architecture is that the capacity of the switches and the number of links between layers increase as one moves up the hierarchy, ensuring that the network can handle high volumes of traffic without bottlenecks.

Routing

Static routing is majorly deployed on fat-tree systems. Static routing forwards packets in a deterministic manner such that the packets between the same source-destination node pair always follow the same routing path. The "destination modulo k" (D-mod-k) algorithm is a widely adopted static routing technique that effectively distributes traffic across links in a fat-tree topology, demonstrating strong performance characteristics [START_REF] Gomez | Deterministic versus adaptive routing in fat-trees[END_REF]. This method determines the next upward link at each hierarchical level by utilizing the destination node's ID, continuing this process until it reaches the shared ancestor node. Subsequently, the algorithm selects the appropriate downward links to direct the traffic towards the intended destination.

In contrast to static routing, adaptive routing forwards packets based on real-time network estimation, thus packets from the same message can be delivered using different paths [START_REF] Kang | Q-adaptive: A multi-agent reinforcement learning based routing on dragonfly network[END_REF][START_REF] Jiang | Indirect adaptive routing on large scale interconnection networks[END_REF]. In addition to static routing, several works also bring adaptive routing to fattree topologies [START_REF] Smith | Mitigating inter-job interference using adaptive flow-aware routing[END_REF][START_REF] Rocher-González | Adaptive routing in infiniband hardware[END_REF]. Adaptive routing allows packets between a source-destination node pair to follow different forwarding paths. While routing flexibility is improved, special designs for choosing the best path and handling out-of-order packet delivery need to be implemented.

Discussions

Fat-tree topology is a popular choice for HPC and large data center networks due to is hierarchical design. However, there are several points that need to be carefully considered for deploying fat-tree. Especially the ones discussed below.

• Scalability: One of the most significant advantages of fattree topology is its scalability. Due to its hierarchical nature and the increased bandwidth at higher layers, fat-trees can easily accommodate growth in network size by adding more nodes at the edge or enhancing the capacity of the upper layers. • High Bandwidth: The increased link capacity and parallel paths available in fat-tree topology allow for high bandwidth between any two nodes in the network. This is particularly beneficial for applications requiring significant data transfer, such as cloud computing and big data analytics. • Fault Tolerance: Fat-tree topology offers inherent fault tolerance through multiple redundant paths between any two points in the network. If a link or switch fails, the network can dynamically reroute traffic, minimizing downtime and ensuring continuous availability. • Complexity: The design and implementation of a fat-tree topology can be complex, especially as the network scales.

Managing the numerous switches and potential paths requires sophisticated network management tools and algorithms to ensure efficient traffic flow and fault management. • Cost: The initial setup cost of a fat-tree topology can be high due to the need for switches with higher bandwidth capacities, especially at the upper layers of the hierarchy.

Operational costs can also increase with the need for advanced network management and monitoring tools.

Fat-tree topology offers a scalable, high-bandwidth, and faulttolerant network architecture that is well-suited to the demands of modern high-performance computing environments. While it presents challenges in terms of complexity and cost, the benefits it provides in terms of performance and reliability make it an attractive choice for data centers, cloud computing infrastructure, and supercomputing.

DRAGONFLY TOPOLOGY

Dragonfly is one of the most popular topologies used in the design of high-performance computing networks and supercomputers. The topology demonstrates its innovation and efficiency by deploying on the two most powerful HPC systems in the Top 500 November list[17].

Architecture

Dragonfly is specifically designed to address the challenges of scalability, low latency, and high bandwidth requirements from nowadays large-scale workloads and applications. The name "dragonfly" is inspired by the natural structure of a dragonfly's wings, which symbolizes the topology's flexibility to connect multiple groups of nodes for different system scales.

As shown in Figure 2, the dragonfly topology is a hierarchical network architecture that consists of groups of routers. These groups are interconnected in a way that minimizes the number of hops data must travel between any two points in the network. This design significantly reduces latency and increases the overall efficiency of data transmission [START_REF] Kim | Technology-driven, highly-scalable dragonfly topology[END_REF]. Each group in a dragonfly network contains a set of routers, and each router is connected to routers within the same group and to routers in other groups using local channels and global channels respectively. This creates a high level of redundancy and paths for data to travel, enhancing the network's fault tolerance and reliability. Many flavors of dragonfly have been deployed depending on their intra-group connection. For example, the current Slingshot network deployed on Frontier and Aurora supercomputer constructs a one-dimensional full intragroup network [START_REF]Aurora supercomputer[END_REF][19] [START_REF] Sensi | An in-depth analysis of the slingshot interconnect[END_REF]. Aries network arranges routers in a group as a two one-dimensional grid, and each dimension is fully connected [START_REF] Alverson | Cray xc series network[END_REF]. Dragonfly+ constructs routers into a tree-like bipartite graph intra-group connection [START_REF] Shpiner | Dragonfly+: Low cost topology for scaling datacenters[END_REF][23].

Advantages of Dragonfly Topology

One of the primary advantages of the dragonfly topology is its scalability. As computational demands grow, networks designed with dragonfly topology can be expanded by adding more groups or increasing the number of routers within groups, without significantly impacting the network's performance. This makes it an ideal choice for supercomputing environments where computational needs are continuously evolving.

Another significant advantage is the high bandwidth and low latency guaranteed by different dragonfly flavors [START_REF] Kang | Modeling and analysis of application interference on dragonfly+[END_REF][START_REF] Salimi | An analysis of long-tailed network latency distribution and background traffic on dragonfly+[END_REF]. By minimizing the number of hops between nodes, the dragonfly topology ensures that data packets reach their destination more quickly compared to traditional network designs. This is particularly beneficial in applications that require real-time data processing and transmission, such as scientific simulations and data analytics.

Challenges and Considerations

While the dragonfly topology offers numerous benefits, it also presents certain challenges. The complexity of the network's design can make it difficult to implement and manage, especially in very large-scale systems. Additionally, the initial cost of setting up a dragonfly network can be higher than traditional networks due to the specialized hardware required. Another challenge in network design that should not be overlooked is the network workload interference [START_REF] Kang | Workload Interference Analysis and Mitigation on Dragonfly Class Networks[END_REF][START_REF] Pu | Understanding performance interference of i/o workload in virtualized cloud environments[END_REF], which usually requires additional QoS or congestion control mechanisms to provide a higher network capability [START_REF] Sensi | An in-depth analysis of the slingshot interconnect[END_REF][START_REF] Jeremiah | Opportunities and limitations of qualityof-service in message passing applications on adaptively routed dragonfly and fat tree networks[END_REF][START_REF] García | Ofar-cm: Efficient dragonfly networks with simple congestion management[END_REF]. However, these challenges are often outweighed by the long-term benefits of improved performance and scalability.

In conclusion, the dragonfly topology represents a significant advancement in the design of network architectures for high-performance computing environments. Its ability to provide low latency, high bandwidth, and exceptional scalability makes it an ideal choice for supercomputers and large-scale data centers. Despite its challenges, the dragonfly topology continues to gain popularity as the demand for faster and more efficient computing networks grows.

HYPERX TOPOLOGY

HyperX topology is a significant concept in the field of network architecture, particularly relevant to the design of high-performance computing (HPC) systems and data centers [START_REF] Ho Ahn | Hyperx: topology, routing, and packaging of efficient large-scale networks[END_REF]. Unlike the more commonly known dragonfly topology, which is characterized by its hierarchical structure, HyperX is recognized for its flexibility and scalability, offering an alternative approach to constructing large-scale, high-bandwidth, and low-latency networks [START_REF] Azizi | On the topological properties of hyperx[END_REF].

Architecture

HyperX topology is designed to be a highly scalable and efficient network architecture that can support a vast number of nodes with relatively low latency and high bandwidth. It is a generalization of the flattened butterfly and flattened torus topologies, allowing for more flexibility in network design and scalability. The HyperX topology constructs high-dimensional interconnection networks by crossing multiple layers of the network, which can be tailored to the specific needs and scale of a computing environment.

The key characteristic of HyperX topology is its multi-dimensional mesh structure, which allows for multiple paths between any two nodes in the network. In HyperX, routers on the same dimension are fully connected. For example, figure3 illustrates the link connection of two HyperX networks at different dimensionality. Figure3(a) shows a two-dimensional HyperX with two routers in the first dimension and four routers in the second. Figure3(b) depicts a three-dimensional HyperX, with three routers on each dimension. The redundancy in paths between two nodes enhances the network's fault tolerance and load balancing capabilities, as traffic can be rerouted in the event of a node or link failure, ensuring consistent performance and reliability.

Advantages of HyperX Topology

One of the primary advantages of the HyperX topology is its scalability. As computational requirements grow, networks designed with HyperX can easily be expanded by adding more dimensions or nodes within each dimension. This flexibility makes it an attractive option for evolving HPC environments and data centers that need to accommodate increasing data volumes and computational demands.

Moreover, HyperX offers high bandwidth and low latency by providing direct paths between nodes and reducing the number of hops required for data transmission [START_REF] Kathareios | Cost-effective diameter-two topologies: Analysis and evaluation[END_REF]. This is particularly beneficial for applications that require rapid communication between nodes, such as parallel computing tasks and real-time data processing.

Challenges and Considerations

Despite its advantages, implementing HyperX topology comes with its set of challenges [START_REF] Domke | The first supercomputer with hyperx topology: A viable alternative to fat-trees?[END_REF][START_REF] Kang | Workload interference prevention with intelligent routing and flexible job placement on dragonfly[END_REF]. Designing and managing a multidimensional network can be complex, requiring advanced planning and specialized knowledge to optimize performance and ensure reliability [START_REF] Mudigonda | Taming the flying cable monster: A topology design and optimization framework for {Data-Center} networks[END_REF]. Additionally, the initial setup and hardware costs can be higher compared to traditional network topologies, making it a significant investment.

In conclusion, HyperX topology represents a powerful and flexible approach to designing network architectures for high-performance computing environments. Its scalability, high bandwidth, and low latency make it an ideal choice for a wide range of applications, from scientific simulations to cloud computing. While there are challenges in implementing HyperX topology, its potential benefits in terms of performance and adaptability make it a compelling option for future network designs in HPC systems and beyond. As computational demands continue to grow, the importance of efficient and scalable network topologies like HyperX will undoubtedly increase, driving further innovation in this field.

CONCLUSION

In this paper, we analyzed three popular interconnect topologies that are widely adopted in HPC and data center large-scale parallel systems. The studied topologies are fat-tree, dragonfly, and HyperX networks. All these three topologies show their advantages in supporting extreme-scale systems with high bandwidth, low latency, and flexibility to support systems at different scale.

Figure 1 :

 1 Figure 1: An example of a three-level fat-tree

Figure 2 :

 2 Figure 2: Hierarchical connections in dragonfly topology. (a) groups are all-to-all connected through the inter-group network with global channels. (b) routers within a group are linked through local channels and connected to the terminals through terminal channels.

Figure 3 :

 3 Figure 3: Link connections between routers in the HyperX topology, nodes are omitted. (a) Two-dimensional HyperX with two and four routers on each dimension. (b) Three-dimensional HyperX, with three routers on each dimension.


Corentin Ségalas
email: corentin.segalas@u-bordeaux.fr

Catherine Helmer

Robin Genuer

Cécile Proust-Lima

Functional principal component analysis as an alternative to mixed-effect models for describing sparse repeated measures in presence of missing data

Keywords: Functional principal component analysis, sparse functional data, missing data, mixed models

Analyzing longitudinal data in health studies is challenging due to sparse and error-prone measurements, strong within-individual correlation, missing data and various trajectory shapes. While mixed-effect models (MM) effectively address these challenges, they remain parametric models and may incur computational costs. In contrast, Functional Principal Component Analysis (FPCA) is a non-parametric approach developed for regular and dense functional data that flexibly describes temporal trajectories at a lower computational cost. This paper presents an empirical simulation study evaluating the behaviour of FPCA with sparse and error-prone repeated measures and its robustness under different missing data schemes in comparison with MM. The results show that FPCA is well-suited in the presence of missing at random data caused by dropout, except in scenarios involving most frequent and systematic dropout. Like MM, FPCA fails under missing not at random mechanism. The FPCA was applied to describe the trajectories of four cognitive functions before clinical dementia and contrast them with those of matched controls in a case-control study nested in a population-based aging cohort. The average cognitive declines of future dementia cases showed a sudden divergence from those of their matched controls with a sharp acceleration 5 to 2.5 years prior to diagnosis.

Introduction

In medical research, it is very common to collect in a cohort of participants the repeated measures of various markers at different visits in order to study the trajectories of the underlying biological processes. For example, in HIV research studies, the temporal evolution of CD4 is used to assess the immune system response [START_REF] Boscardin | Longitudinal models for AIDS marker data[END_REF]. In prostate cancer, the prostate specific antigen is a marker used to monitor prostate cancer progression [START_REF] Ilic | Prostate cancer screening with prostate-specific antigen (PSA) test: a systematic review and meta-analysis[END_REF]. In cognitive ageing studies, the temporal evolution of neuropsychological tests is used to assess the cognitive abilities trend among the older adults [START_REF] Alpérovitch | Les études épidémiologiques sur le vieillissement en France : de l'étude Paquid à l'étude des Trois Cités[END_REF].

Marker data repeatedly measured over time in longitudinal studies have several characteristics that raise analytical challenges. First, the structure of the data induces a strong within-individual correlation which implies a variability across individuals. Second, the data collection process [START_REF] Lokku | Summarizing the extent of visit irregularity in longitudinal data[END_REF], often heavy in medical studies, usually leads to sparse measurements and individual-specific timings. For instance, in epidemiological cohorts of chronic diseases, it is not rare to have data every 2 to 5 years only [START_REF] Boscardin | Longitudinal models for AIDS marker data[END_REF][START_REF] Ilic | Prostate cancer screening with prostate-specific antigen (PSA) test: a systematic review and meta-analysis[END_REF][START_REF] Alpérovitch | Les études épidémiologiques sur le vieillissement en France : de l'étude Paquid à l'étude des Trois Cités[END_REF]. Third, participants may leave the study before the end because they died or because they refused to undergo the next visit which leads to dropout [START_REF] Hogan | Handling drop-out in longitudinal studies[END_REF]. These early truncations of the longitudinal process are often linked to the health status of the participant, and thus are likely informative. Fourth, marker data collected in health studies are systematically noisy observations of the true underlying process of actual interest. Finally, the marker trajectories may exhibit different shapes, requiring flexible modeling tools of time trends. Figure 1 illustrates the characteristics of such marker data in the context of four neuropsychological tests collected in a cognitive aging study with the error-prone trajectories and sparse measurements (left panels), the individual-specific timings (middle panels) and the missingness patterns (right panels).

Mixed-effect models (MM) are the standard method to analyze longitudinal data in health studies [START_REF] Verbeke | Linear Mixed Models for Longitudinal Data[END_REF]. They extend the classic lineal regression models to correlated data by modeling both marginal trajectories using fixed effects and the individual-specific deviations using individual regression coefficients called the random effects. MM split the observations into the observation error and the underlying data generation process defined in continuous time, thus handling any measurement time and accounting for serial correlation. Any basis of time functions can be considered, including splines or polynomials, for flexible modeling. As estimated in the Maximum Likelihood framework, MM estimates are robust to missing data as long as the missingness mechanism is predictable by the observed data (mechanism also known as Missing at Random) [START_REF] Rubin | Inference and Missing Data[END_REF][START_REF] Roderick | Statistical Analysis With Missing Data[END_REF]. MM thus handle most challenges faced with the longitudinal data collected in health data. However, MM are parametric tools; the distribution of both the error and the random effects must be specified by the user. Also, in some specific cases, the flexibility put into the MM may come with a great computational cost, particularly when the model is nonlinear and includes many random effects because of an heavy numerical integration in the log-likelihood calculation [START_REF] Davidian | Nonlinear models for repeated measurement data: An overview and update[END_REF].

Repeated measures of markers are observations of a continuous-time process and as such they could also be considered as sparse and irregular functional data. Indeed, classic functional data is a collection of noisy observations from an unknown random function f on a regular and temporally dense grid of time. In functional data analysis, the main goal is to describe and summarize the temporal evolution of these trajectories. One of the most well-known method to do so is the functional principal component analysis (FPCA) [START_REF] Ramsay | Functional Data Analysis[END_REF][START_REF] Wang | Review of Functional Data Analysis[END_REF]. Similarly to classic principal component analysis, it projects functional data into a lower dimensional functional space defined by a functional basis. In this functional space, each participant has its own set of coordinates called functional principal component scores. In health studies, FPCA is also very appealing because it allows, with a fully non parametric approach, to describe and summarize temporal trajectories in a parsimonious way using only a few principal functional components and individual contributions to these components. However, to be suitable for the analysis of longitudinal data stemmed from health studies, functional analysis methods should necessarily handle the challenges of sparsity, irregularity of the observation grid and missing data. Previous work have proposed FPCA for sparse longitudinal data [START_REF] Yao | Functional Data Analysis for Sparse Longitudinal Data[END_REF]. Furthermore, FPCA should work in theory with irregular observation grids although current implementations in R packages are not always flexible enough. The main unsolved challenge is the behaviour of FPCA techniques with missingness mechanisms. No theoretical or numerical evaluation has been sought. While FCPA is expected to behave well in case of completely random missingness since it constitutes a specific type of irregular grid, the behaviour of the FPCA to other types of missingness still deserve to be investigated, especially when working with longitudinal data from health studies where informative missingness is common. The aim of this paper is to evaluate the robustness of FPCA to various schemes of missing data in an empirical simulation study and to compare its performances to the classic mixed-effect model approach.

Section 2 briefly introduces the motivating data example with repeated evaluations of cognitive functioning in a population-based aging cohort. In Section 3, the statistical framework -in particular the missing data mechanisms -is described and both the mixed-effect models and the FPCA are introduced. Results from two empirical simulation studies are reported in Section 4. The two methods are then applied in Section 5 to describe the trajectories of four cognitive functions in years preceding clinical dementia onset in contrast with the trajectories of matched controls. Our findings are summarized and discussed in Section 6.

The Three-City Study

The Three-City Study (3C Study) is a French observational cohort study started in 1999 aiming at exploring the link between vascular diseases and dementia in the older adults [START_REF] Alpérovitch | Les études épidémiologiques sur le vieillissement en France : de l'étude Paquid à l'étude des Trois Cités[END_REF]. Subjects aged 65 years and older at baseline were recruited in three French cities and followed-up every 2 to 3 years during 17 years. At each visit, the cognitive functioning was measured by a battery of cognitive tests, and a diagnosis of dementia was established according to a three steps procedure: a screening based on neuropsychological performance; a clinical evaluation performed by a neurologist; and a final diagnosis made by an independent committee of neurologists (see [START_REF] Alpérovitch | Les études épidémiologiques sur le vieillissement en France : de l'étude Paquid à l'étude des Trois Cités[END_REF] for details).

Our motivating aim was to describe the cognitive trajectories of participants diagnosed with dementia and contrast them with those expected during natural aging. To this end, we built a nested case-control study from the 3C Study. We included the 174 incident dementia cases diagnosed in the Bordeaux Center with at least one measure at each cognitive test, and individually matched them in density of incidence with 156 controls free of dementia at the time of diagnosis of the case. The matching was done according to sex, presence of at least one APOE-ϵ4 allele (main genetic risk factor of dementia), level of education, age at diagnosis with a 3-year margin, and length of follow-up at diagnosis (in years). We removed the cognitive measures occurring 5 years after the diagnosis (considered as time 0) as most cases were not able anymore to undergo a large cognitive battery.

We focused on four psychometric scores that evaluate different cognitive aspects: the Mini-Mental State Examination score (MMSE) [START_REF] Folstein | Mini-mental state". A practical method for grading the cognitive state of patients for the clinician[END_REF] which evaluates global cognitive functioning (with range 0 to 30), the Benton Visual Retention Test score (BVRT) [START_REF] Benton | The revised visual retention test: clinical and experimental applications. The Psychological Corporation[END_REF] which measures visual memory (with range 0 to 15), the 30-second Isaacs Set Test (IST) score [START_REF] Isaacs | The Set test as an aid to the detection of dementia in old people[END_REF] which assesses verbal fluency (with range 0 to 73) and the number of correct moves per minute at the Trail Making Test score part A (TMTA) [START_REF] Reitan | Validity of the Trail Making Test as an Indicator of Organic Brain Damage[END_REF] which assesses executive functioning and attention. Note that the TMTA was not assessed at the second visit.

Figure 1 describes the sample with the individual cognitive trajectories, the timings and patterns of missingness, thus illustrating the sparsity, irregularity in timings and missing data usually encountered in cohort studies.

Methods

Sparse and irregular repeated measures data

We assume that we observe repeated measures of a longitudinal marker over time t for each unit i of a sample of size N . These observations y ij also noted y i (t ij ) with j = 1, . . . , n i lead to an individual vector of observations y i = (y i1 , . . . , y ini ) ⊤ where n i denotes the number of observations for unit i. The observations y ij are measurements with error of the underlying process of interest y ⋆ i (t) so that y ij = y ⋆ (t ij ) + ε ij where ε ij denotes the independent random observation error.

Classically, in the functional data framework, the observation grid t = (t 1 , . . . , t n ) ⊤ is dense, regular over time (e.g. hourly measures, daily measures) and the same across the n units. In contrast, in health studies on humans (cohorts, clinical trials), the data collection is frequently done with a lower frequency and the visit times t i = (t i1 , . . . , t ini ) ⊤ may vary across individuals. For example, in the 3C cohort data, there may be variations of several months or even years across individuals around the theoretical visits every 2 to 3 years because collecting cognitive test results implies an heavy and expensive observation process. Additionally, it is very common to have missing data which leads to different size of observation grid. Such a design defines sparse and irregular functional data.

Missing data mechanisms

We use the classical formalism of missing data for a marker of interest as introduced by Little and [START_REF] Roderick | Statistical Analysis With Missing Data[END_REF] [START_REF] Roderick | Statistical Analysis With Missing Data[END_REF]. We denote Y the set of all marker data which includes those that are observed Y o and those that are missing Y m . The binary variable R is the indicator variable of observation with R ij = 1 if Y ij is observed, and 0 otherwise. Little and [START_REF] Roderick | Statistical Analysis With Missing Data[END_REF] [START_REF] Roderick | Statistical Analysis With Missing Data[END_REF] distinguishes three mechanisms of missingness. Data is called Missing Completely At Random (MCAR) when the probability that the marker is missing does not depend on the marker at all,

P (R|Y o , Y m ) = P (R).
It happens for instance when a cognitive test is unavailable at a certain visit because this test was not scheduled in the protocol for this visit (as this is the case for the TMTA at the second visit in 3C). Data is called Missing At Random (MAR) when the probability that the marker is missing depends on the observed values of the marker,

P (R|Y o , Y m ) = P (R|Y o ).
It happens for instance when a participant refusal to pass a cognitive test depends upon his/her results at previous visits. Finally, data is called Missing Not At Random (MNAR) when it can also depend on the unobserved values of the marker,

P (R|Y o , Y m ) = P (R|Y o , Y m ).
This can happen for instance when a participant refuses to pass the test because he/she feels that his/her cognitive level has dropped since last visit. Missing data can be divided into two main types: intermittent missing data, when a measurement is missing during the follow-up but other measurement can be collected later on, and dropout data which happens when any visit after the dropout visit will be missing. In this work, we are interested in both types of missing data.

Mixed-effect models

Mixed-effect models have been introduced by Laird and Ware (1982) [START_REF] Laird | Random-Effects Models for Longitudinal Data[END_REF] and since then have been widely used to model repeated measures collected in longitudinal health studies. We focus on a mixed model that assumes a continuous marker with zero-mean independent Gaussian errors ε ij so that

y ij = y ⋆ i (t ij ) + ε ij .
The underlying value of the marker at time t is then modelled as follows:

y ⋆ i (t) = g(t, X i (t), β, b i ) (1)
where X i (t) denotes a set of time-dependent covariates, β denotes the vector of fixed parameters, b i denotes the vector of random effects with b i ∼ N (0, B) and g a basis of functions that describes the shape of the marker trajectory. The observation error ε ij ∼ N (0, σ 2 ) is assumed to be independent from the vector of random effects. In the reminder of this work, we concentrate on linear mixed-effect models where

g(t, X i (t), β, b i ) = p k=1 (β k +b ik )g k (t, X i (t)
) with (g 1 , . . . , g p ) a p-basis of a priori specified functions of time, most often splines or polynomials.

Mixed-effect models have been extended to handle different types of data: binary or count data using generalized linear mixed model theory, or continuous data with strong ceiling or floor effect which is very common among neuropsychological scores [START_REF] Proust-Lima | Sensitivity of four psychometric tests to measure cognitive changes in brain aging-population-based studies[END_REF][START_REF] Proust-Lima | Misuse of the linear mixed model when evaluating risk factors of cognitive decline[END_REF]. For the latter, a parametric monotonic function is applied to Y in Equation 1 to normalize the marker data.

The mixed-effect model is generally estimated by maximizing the log-likelihood:

ℓ(θ) = N i=1 log f (Y i |b i )f (b i )db i .
where θ denotes the model parameters (θ = (β, vec(B), σ) ⊤ ), f (Y i |b i ) is the density function of the outcome conditional to the random effects and f (b i ) is the density of the random effects. When the model is linear in the random effects, the integral has an analytical solution. In other cases, the loglikelihood computation often implies to numerically approximate the integral [START_REF] Pan | Gauss-Hermite Quadrature Approximation for Estimation in Generalised Linear Mixed Models[END_REF]. Fitted values of existing observations or predictions of new data are easily derived from a fitted mixed model using the a posteriori distribution of the random effects given the individual observations. In linear mixed-effect models, it is directly computed using conditional expectation properties of multivariate normal variables while for nonlinear models, it is approximated by the mode of the posterior distribution [START_REF] Laird | Random-Effects Models for Longitudinal Data[END_REF][START_REF] Molenberghs | Linear Mixed Models for Longitudinal Data[END_REF]. Despite the possible use of flexible families of time-functions (e.g., splines), a MM remains a parametric model with a distribution for the random effects and for the random error to be chosen (often Gaussian). Additionally, some assumptions on the variance structure of the random effect might be needed when the number of random effects increases to avoid a too large number of covariance parameters.

One essential feature of the MM is its robustness to MAR data, which ensures an unbiased estimator under MAR [START_REF] Rubin | Inference and Missing Data[END_REF]. When the missing data are suspected to be MNAR, that is when the missing mechanism also depends on unobserved data, a joint model of the missingness mechanism is required. This can be done with a logistic model for intermittent missing data or a survival model for dropout [START_REF] Diggle | Informative Drop-Out in Longitudinal Data Analysis[END_REF][START_REF] Rizopoulos | Joint Models for Longitudinal and Time-to-Event Data: With Applications in R[END_REF]. However, in any case, the model requires a correct specification of the missing data mechanism to yield robust estimates [START_REF] Thomadakis | Longitudinal and time-to-dropout joint models can lead to seriously biased estimates when the drop-out mechanism is at random[END_REF]. For instance, if the instantaneous risk of dropout is assumed to be associated with the marker level, the following proportional hazard model can be jointly estimated with the mixed model in [START_REF] Alpérovitch | Les études épidémiologiques sur le vieillissement en France : de l'étude Paquid à l'étude des Trois Cités[END_REF]:

α i (t) = α 0 (t) exp( X⊤ ij γ + y ⋆ i (t) ⊤ η)
with α 0 (t) the baseline risk function, X⊤ ij a set of covariates, γ and η two vectors of parameters to be estimated.

Functional Principal Component Analysis

Functional Principal Component Analysis has been introduced by Besse and Ramsay (1986) [START_REF] Besse | Principal components analysis of sampled functions[END_REF] and since then has been used for describing, modeling, predicting or classifying abundant temporal series [START_REF] Lin | A survey of functional principal component analysis[END_REF][START_REF] Wang | Review of Functional Data Analysis[END_REF]. As for mixed-effect models, we consider additive random errors such that observations

y ij = y ⋆ i (t ij ) + ε ij with y ⋆ i (t ij
) the underlying process of interest. In functional data, the vector y ⋆ i = {y ⋆ (t ij )} j=1,...,ni for unit i is considered as a random realization of an unknown underlying function. Consequently, the whole set of longitudinal trajectories y ⋆ = (y ⋆ 1 , . . . , y ⋆ N ) ⊤ is a collection of realizations of an unknown random function f assumed to be smooth. This unknown function has mean µ(t) and covariance function G(s, t) = cov(f (t), f (s)). FPCA is essentially based on the Karhunen-Loève decomposition [START_REF] Karhunen | Zur Spektraltheorie stochastischer Prozesse[END_REF][START_REF] Loève | Fonctions aléatoires à décomposition orthogonale exponentielle[END_REF] which states that the underlying unknown smooth function f can be decomposed into:

f i (t) = µ(t) + ∞ k=1 ξ ik ϕ k (t) (2)
where the index i in f i emphasizes the fact that the random realisation of f is subject specific, ξ ik and ϕ k are the eigenvalues (or functional principal component scores) and eigenfunctions from the decomposition of the covariance operator G, respectively. From this decomposition, dimension reduction is straightforward by limiting the sum in (2) to the first K components only:

f i (t) = µ(t) + K k=1 ξ ik ϕ k (t). (3)
From the noisy data y ij , local linear smoothers [START_REF] Ramsay | Functional Data Analysis[END_REF][START_REF] Yao | Functional Data Analysis for Sparse Longitudinal Data[END_REF] are used to estimate the mean function µ and the covariance function G. The estimated covariance surface is discretized to make spectral decomposition of the covariance operator easier and from that, eigenfunctions are obtained [START_REF] Rice | Estimating the Mean and Covariance Structure Nonparametrically When the Data are Curves[END_REF][START_REF] Wang | Review of Functional Data Analysis[END_REF]. In the case of dense functional data, the scores can be estimated using a numerical integration [START_REF] Ramsay | Functional Data Analysis[END_REF] which is not feasible in the case of sparse data. Functional Principal Components Analysis through Conditional Expectation (PACE) [START_REF] Yao | Functional Data Analysis for Sparse Longitudinal Data[END_REF] has been proposed as an alternative way of computing the eigenvalues in case of sparse data. Using properties of conditional Gaussian multivariate density, an expression for E(ξ i |Y i ) is obtained, very close in spirit to what is done in linear mixed-effect models. From that, by plugging in the estimated eigenvalues, eigenfunctions, covariance and mean function, an estimate of the score can be computed. In PACE, the scores are assumed to be Gaussian but the method is robust to violation of this assumption [START_REF] Yao | Functional Data Analysis for Sparse Longitudinal Data[END_REF].

When using FPCA, the determination of K in (3) is critical. Apart from rare exceptions, K is determined from the data by optimizing a statistical criterion such as percentage of variance explained, AIC, or BIC where cross-validation can be used. Fitted values of existing observations or predictions of new data can be easily computed by plugging into (3) the estimated μ(t ij ) and φk (t ij ), and the individual eigenvalues ξik computed using the PACE algorithm.

Many R packages have implemented FPCA for functional data but only a few are suited to the case of sparse and irregular functional data, in particular during the prediction step for new units. For example, fda [START_REF] Ramsay | fda: Functional Data Analysis[END_REF] does not implement PACE and requires the same observation grid for all units which is not realistic in many health contexts. Package face [START_REF] Xiao | face: Fast Covariance Estimation for Sparse Functional Data[END_REF] is suited to sparse and irregular functional data but does not accept missing values in the input data and therefore cannot predict missing observations. Package MFPCA [START_REF] Happ | Multivariate functional principal component analysis for data observed on different (dimensional) domains[END_REF] is suited to sparse and irregular functional data but it struggles to make predictions on new units with different observation grids. In our experience, fdapace [START_REF] Zhou | fdapace: Functional Data Analysis and Empirical Dynamics[END_REF] is the only R package that implements PACE and allows flexible prediction for new units, and was thus used in this work.

Contrarily to linear mixed-effect models, no theoretical justification exists regarding FPCA robustness to missing at random data, and no empirical study has explored this issue so far.

Simulation study

We numerically evaluated the performances of FPCA to analyze sparse, noisy and irregular repeated data in the presence of dropout with (i) a comparison of the predicted values with those obtained by mixed models considered as the gold standard approach, and (ii) an evaluation of the robustness of FCPA estimations to dropout. We followed the recommendations of [START_REF] Morris | Using simulation studies to evaluate statistical methods[END_REF] [START_REF] Morris | Using simulation studies to evaluate statistical methods[END_REF] to plan and report these two simulation studies.

Comparison between FPCA and MM

Aims

The first simulation aimed to contrast the performances of MM and FPCA to predict missing values under different scenarios of dropout.

Data generation mechanism

We generated individual repeated data in samples of N = 700 participants. The visit process was generated from a fixed grid of visits between t = 0 and t = 12 with either 5, 7 or 13 measures per participant (i.e. measures every 3, 2 or 1 year, respectively). To mimic a more realistic observation window, the final visit times t ij were obtained by adding a random uniform noise (within ± 1.5, 1, 0.5 years respectively to measure frequency) around the theoretical visits for each observation. The repeated marker data at each visit time t ij were generated according to a mixed model with a nonlinear shape over time captured by a five-parameter logistic model as implemented in the nlraa package [START_REF] Miguez | nlraa: Nonlinear Regression for Agricultural Applications[END_REF]:

y ⋆ i (t ij ) = y ∞i + y 0i -y ∞i (1 + exp(α i × (log t ij -log τ i ))) θi
where y 0i = 10 + b i1 and y ∞i = 35 + b i2 are the individual-specific asymptotic values of y for 0 and infinite times, respectively; τ i = 10 + b i3 is the time when the mean value of y is attained (only if θ i = 1), i.e. y(τ i ) = (y 0i + y ∞i )/2, α i = 5 + b i4 controls the steepness of the transition and θ i = 1 + b i4 /3 is an asymmetry parameter. Individual-specific trajectories were generated assuming random-effects

b i = (b 1i , b 2i , b 3i , b 4i ) ⊤ ∼ N (0, B) with B a diagonal 4 × 4-matrix with diagonal (3, 5, 2, 2) ⊤ .
The proneto-error observations were derived by adding an independent homoscedastic Gaussian error ε ij ∼ N (0, 3 2 ) such that: y(t ij ) = y ⋆ i (t ij ) + ε ij . At this stage, we have a complete set of repeated measures for each participant up to the administrative censoring at year 12. A right-truncation of the longitudinal process was added by considering a dropout (i.e., early leave from the study) with two intensities (30% and 60%) and the six following missing data scenarios (see supplementary material for details):

• MCAR: the participant drops out at t ij with a probability determined by a logistic model with t ij as predictor; • threshold fixed MAR: the participant systematically drops out at t ij if Y i(j-1) is above a threshold ν;

• threshold MAR: the participant drops out at t ij with a probability determined by a logistic model with the indicator that Y ij-1 > ν as predictor;

• increasing MAR: the participant drops out at t ij with a probability determined by a logistic model with Y ij-1 as predictor;

• threshold MNAR: the participant drops out at t ij with a probability determined by a logistic model with the indicator that Y ij > ν as predictor;

• increasing MNAR: the participant drops out at t ij with a probability determined by a logistic model with Y ij as predictor;

The generated longitudinal data of 100 randomly selected individuals are displayed in Figure 2. The right panels illustrate the two types of early truncation of the longitudinal process: the dropout under the above mechanisms (with a rate of 0.6) and the administrative censoring at year 12.

Estimands

Each sample was separated into a training set of 200 participants and a test set of 500 participants. The estimands of interest were the predictions of the marker values y ij in the test set whether they were observed or missing (true generated values were known from the generation procedure). These predictions were computed from the estimations of FPCA and MM on the training set.

Methods

For each simulation setting, that is for a fixed frequency of visits (every 3, 2 and 1 year), for a fixed rate of dropout (30% or 60%), for a fixed missing data scenario (MCAR, threshold fixed MAR, threshold MAR, increasing MAR, threshold MNAR, increasing MNAR), we estimated the following models:

• Linear mixed-effect models (including a random effect on the intercept and on each time function)

with time trend: modelled as a quadratic trend (LMM quad); modelled as a cubic trend (LMM cub); approximated by natural splines with 2 internal knots at the terciles (LMM spl quant); approximated by natural splines with 2 equidistant internal knots (LMM spl equi).

• Functional Principal Component Analysis with the number of principal components: fixed at 4 (FPCA k4); chosen such that the percentage of explained variance reached 90% (FPCA fve90) or 99% (FPCA fve99).

• Shared random effect joint model (JM) in which the longitudinal submodel followed a quadratic time trend (random intercept and random slope on each time function) and the survival submodel was a proportional hazard model with a baseline risk function approximated by B-splines (3 internal knots at the quartiles) and a dependence through the current underlying value of the marker. Note that the joint model was only estimated in the two MNAR scenarios (threshold MNAR and increasing MNAR) as a gold standard comparator for the MNAR situation.

The linear mixed-effect models were estimated using the lcmm package [START_REF] Proust-Lima | lcmm: Extended Mixed Models Using Latent Classes and Latent Processes[END_REF] (function hlme) and predictions were obtained using the predictRE function. The joint models were estimated using the JM package [START_REF] Jm | An R package for the joint modelling of longitudinal and time-to-event data[END_REF], and predictions computed with the predict function. FPCA models were estimated with the fdapace package [START_REF] Zhou | fdapace: Functional Data Analysis and Empirical Dynamics[END_REF]. In contrast with previous methods, prediction was only possible for times available in the training set. We thus discretized the time window (with a thin step of 0.5) to ensure that every discretized t ij in the test set also appeared in the training set, and used the internal predict function.

Performance measure

For each setting, the procedure was replicated on 1000 replicates. The performances were evaluated with the root mean square error (RMSE) by separating: (i) the RMSE computed on the actual observations of the set test (i.e., non-missing due to dropout), and (ii) the RMSE computed on the missing data of the test set (after time of dropout).

To make all the RMSEs comparable across models within a fixed scenario, we considered relative RMSEs. RMSEs were standardized by the RMSE computed on the same set test according to a reference model estimated in the complete training set (i.e. in the absence of dropout). The reference model was a mixed model with B-splines time trend (3 internal knots placed at the quartiles).

Results

The relative RMSEs of each model computed on the missing test data are reported in Figure 3 (see Supplementary Figures 2, 3 and4 for the observed test data) for four missing data scenarios (MCAR, fixed MAR, increasing MAR and increasing MNAR), two rates of dropout (30% and 60%), and the three frequencies of visits.

As expected, virtually the same relative RMSE were obtained for MM and FPCA in the MCAR case, whatever the frequency of visit and the rate of dropout. With a systematic MAR dropout after an observation above a given threshold, the MM relying on splines showed the lowest relative RMSEs. With a 30% rate of dropout, those of FPCA were comparable or even lower when data were denser (every year). However, with a 60% rate of dropout, the performances of FPCA dropped, probably due to the sudden and local loss of information. With a MAR dropout that increases as the last observed marker value increased, the dropout is more spread along the observation window and FPCA lead to very competing relative RMSEs compared to MM with splines time trends. This is true whatever the rate of dropout, and even provide slightly better fit for denser data and lower rate of dropout, underlying their Relative RMSEs computed on missing test data over 1000 replicates for four missing data scenarios (MCAR, fixed MAR, increasing MAR and increasing MNAR from left to right), two rates of dropout (30% and 60%) and three frequencies of visits (every 3, 2 and 1 year from top to bottom). Each boxplot summarises the RMSEs of a model relative to the reference model (from left to right): LMM quad, LMM cub, LMM spl quant, LMM spl equi, FPCA fve90, FPCA fve99, FPCA k4, JM (only for the MNAR case).

flexibility to capture complex trends. With a MNAR dropout that increases as the simultaneous marker level increases, MM and FPCA models show larger relative RMSEs compared to the more appropriate JM. However, again FPCA shows worse behaviour than MM when the rate of dropout is high (60%). In summary, this simulation study suggests that FPCA behave relatively similarly as the flexible linear mixed-effect models in MCAR, MAR and MNAR scenarios when the dropout rate remains moderate (30%). In some specific situations with a high rate of dropout, the FPCA seems to miss the underlying pattern of trajectory. These results suggest that FPCA may constitute a suitable alternative to MM in non pathological and not too prevalent MAR settings.

The results for the threshold MAR and threshold MNAR cases, displayed in Supplementary Figure 1, lead to the same conclusions.

Robustness of FPCA to dropout

Aims

This second simulation study aimed to evaluate to which extent the components and scores of FPCA were correctly estimated in the presence of dropout.

Data generation mechanism

The data generation mechanism was exactly the same as in the first simulation study except that the samples were of size N = 200 and the time function for the underlying marker trajectory was generated according to a two-component FPCA rather than a mixed model.

The generating mean function and principal components, displayed in the bottom-left panel of Figure 2, were chosen as those obtained after estimating a FPCA on a dataset simulated as in Section 4.1. For each individual, we then randomly sampled Gaussian scores ξ i ∼ N (0, Ξ), (with Ξ taken as the empirical variance of the scores derived from the FPCA) and applied equation ( 3) to generate the functional data at each observed visit time.

Estimands

The quantities of interest were the estimated mean function and the estimated functional principal components.

Methods

For each setting (rate of dropout, mechanism of dropout and frequency of visits), we applied a FPCA implemented in fdapace package [START_REF] Zhou | fdapace: Functional Data Analysis and Empirical Dynamics[END_REF] to the observed data, fixing the number of components to 2.

Performance measure

We evaluated the relative bias between the estimated mean function or principal components, and their true counterparts on a fixed time grid over 1000 replicates. Additionally, we visually compared the empirical distribution the estimated mean function and functional principal components to the generated ones.

Results

The simulation results are reported in Figure 4 for both the MCAR and increasing MAR scenarios when considering a measure every 2 years. As expected, in the MCAR scenarios (upper panels), there is no bias in the estimations whatever the dropout rate (0%, 30% or 60%) and the estimated curves are close to the true generated mean and principal components. In the increasing MAR (lower panels) and the threshold MAR (Supplementary Figure 12) scenarios, the principal components are overall estimated without bias whatever the dropout rate (0%, 30% or 60%) and the variability of the estimates does not substantially change with the dropout rate. The mean function in contrast shows a small bias at the end of the observation window that increases as the rate of dropout increases. Indeed the mechanism of dropout imposes a larger proportion of missing data by the end of the observation window. These observations are further illustrated by the 10 randomly selected estimated curves which highlight the departure from the true mean function in the 30% and 60% dropout cases. With a systematic MAR dropout after an observation above a threshold, the FPCA performed poorly at the end of the observation window as expected due to the sudden and systematic loss of information. In the increasing and threshold MNAR scenarios (Supplementary Figures 13 and14), the relative bias for the mean function is overall larger and especially at the end of the follow-up. The same results (Supplementary Figures 5-10 and 15-20) were observed for the other frequencies of observations.

Application to cognitive trajectories in prodromal dementia

We applied FPCA and compare it to MM to flexibly describe the trajectories of the four cognitive functions (global functioning, visual memory, verbal fluency, executive functioning) before clinical dementia, and contrast them with those of matched controls in the case-control sample nested within the 3C study. Each cognitive function was measured by a specific score: MMSE for global functioning, BVRT for visual memory, IST for verbal fluency and TMTA for executive functioning. The main characteristics of the nested case-control study sample are reported in Table 1 and individual trajectories and missingness mechanisms are plotted in Figure 1.

The analytical sample comprised 74% of women in both groups; 61% had at least finished secondary school and 25.3% were carriers of at least one APOE-ϵ4 allele. The median age at diagnosis was 84.8 (interquartile range (IQR) = 81.6 -88.7) for cases. The matched index age for controls was 84.8 (IQR = 81.6 -88.8). The four cognitive scores distribution implies a worst cognitive condition for cases with a median systematically lower and an interquartile range extending systematically below those of controls. The number of repeated measures per participant was similar between cases and controls and also among cognitive scores except for TMTA, which had fewer measures.

A separate analysis was conducted among cases and controls. For each subset and each cognitive score we estimated: a FPCA with a targeted percentage of explained variance of 99%, and a linear mixed-effect model in which the cognitive trajectory was modeled with natural cubic splines with two internal knots at the quartiles, and boundary knots at -10 and 2 years to reduce the influence of observations beyond this range. With FPCA as implemented in fdapace, no control was possible for reducing the influence of time periods during which the information is sparser. For both methodologies, 95% confidence intervals were computed: directly from the lcmm output based on the inverse of the negative Hessian matrix for the mixed models, and computing bootstrap confidence intervals with 500 bootstrap samples with the function GetMeanCI for the FPCA.

The estimated marginal trajectories of the cognitive markers obtained by the two methods in each case/control group are displayed in Figure 5. Overall, the estimated mean trajectories and 95% confidence intervals were very similar across methods except for slight differences observed at the end of the observation window. The predicted mean trajectories obtained from the linear mixed-effect model were however much smoother than those estimated by the FPCA. The mean trajectories of the four cognitive functions (visual memory, verbal fluency, global functioning and executive functioning) showed very different trends for controls and cases. Controls had a slow linear decline over time for all the cognitive functions. Far from dementia diagnosis, cases also had a slight linear decline, almost overlapped with the one of the controls for global functioning, very close for visual memory, and almost parallel but approximately 0.1 points and 5 points below for executive functioning and verbal fluency respectively. Then closer to diagnosis, each cognitive decline showed a sharp acceleration among cases: around 5 years before diagnosis for verbal fluency, 3 years before diagnosis for executive functioning and 2.5 years before diagnosis for visual memory and global functioning.

Discussion

Initially developed for dense and regular functional data, FPCA has been recently extended to handle sparsity [START_REF] Yao | Functional Data Analysis for Sparse Longitudinal Data[END_REF]. However, the behaviour of FPCA in the presence of missing data remained uncertain. In this manuscript, we have conducted an empirical simulation study to evaluate the robustness of FPCA in the presence of missing data caused by dropout. Our findings demonstrate that, in practical applications, FPCA is well-suited for analyzing longitudinal data, even when missing data follow a Missing at Random (MAR) mechanism, exhibiting comparable performance to linear mixed-effect models known for their robustness in such scenarios. This observation holds true across various MAR scenarios, observation grid densities, and dropout rates with exception for specific scenarios with more frequent and systematic dropout. In contrast, FPCA along with linear mixed-effect models do not seem to be robust to Missing Not at Random data.

This empirical evaluation specifically targeted dropout, a particular category of missing data pattern wherein the absence of data at any point implies that all subsequent visits are also missing. The dropout scenario may induce a larger information loss than intermittent missing data where information from the longitudinal trajectory can be retrieved, even at the end of the window time. Hence, we anticipate the robustness to dropout found in this work to also apply to intermittent missing data. Note that it is important to distinguish missing data due to dropout and missing data truncated by death [START_REF] Kurland | Longitudinal Data with Follow-up Truncated by Death: Match the Analysis Method to Research Aims[END_REF][START_REF] Rouanet | Interpretation of mixed models and marginal models with cohort attrition due to death and drop-out[END_REF]. In this work, we only focused on missing data due to dropout.

In our application study, without any assumption on the shape of the cognitive trajectories, both methods (based on MM and FPCA) showed very similar trends. Despite the absence of control on the FPCA degree of smoothing using fdapace [START_REF] Zhou | fdapace: Functional Data Analysis and Empirical Dynamics[END_REF], the mean predicted trajectories were smooth over time except for the TMTA for which variations at the very beginning of the observation window suggested an overfitting (see Figure 5). Compared to other scores, TMT was not collected at the second visit of the cohort and fewer measures were available between -15 and -10 (see last row in Figure 1).

Cognitive trajectories exhibited a linear mean trend for controls throughout the window of interest. In contrast, the cognitive trajectories of cases, closely resembling those of controls distant from the diagnosis, exhibited a break with an accelerated decline as the diagnosis approached. This two-stage trend in pathological cognitive decline leading toward dementia had already been identified in cognitive ageing research [START_REF] Jacqmin-Gadda | Random Changepoint Model for Joint Modeling of Cognitive Decline and Dementia[END_REF][START_REF] Amieva | Compensatory mechanisms in higher-educated subjects with Alzheimer's disease: a study of 20 years of cognitive decline[END_REF]. This empirical study highlights the usefulness of FPCA to closely describe, summarize and predict longitudinal data. As a non-parametric approach, it avoids dependence upon strong assumptions and can adapt to any shape of trajectory. However, inference tools are not straightforwardly available with FPCA contrarily to MM which benefits from the maximum likelihood framework. For instance, parametric changepoint mixed models have been proposed to conduct inferential tasks to specifically explore this cognitive decline two-stage trend [START_REF] Dominicus | A random change point model for assessing variability in repeated measures of cognitive function[END_REF][START_REF] Segalas | A hypothesis testing procedure for random changepoint mixed models[END_REF].

A parallel between the mixed-effect model and the FPCA frameworks exists. The role of the individual random effects in the mixed-effect model is very similar to the role of the individual scores in the FPCA. They both represent individual deviations to the mean trajectory: the marginal model in mixed-effect model and the mean function in the FPCA. Moreover, under Gaussian assumption, the random effect prediction for linear mixed-effect model is based on Gaussian conditional expectation properties and its formula is very close to FPCA score's formula from the PACE algorithm. Despite this parallel, FPCA and more broadly any method from functional data analysis were initially developed for dense and regular functional data. This work opens up possibilities regarding the utilization of functional data analytical tools in the context of sparsely and irregularly measured health data. However, their application should be conducted carefully and their evaluation remains an active area of research.

Figure 1 :

 1 Figure1: Description of the repeated data in 174 dementia cases and 174 matched controls (N=348) from the 3C study according to the delay (in years) to dementia diagnosis for the cases, and the delay to the dementia diagnosis of the matched case for controls (with 0 corresponding to the dementia diagnosis visit of the case). Are reported the individual observed trajectories over time to diagnosis of four cognitive scores (MMSE, BVRT, IST30 and TMT-A) with 10 randomly selected individual trajectories emphasized in green (first column); the distribution of the times of measurement (second column) and the missingness pattern vizualisation for each individual at each visit (third column). A black segment means the score is observed, a white segment means no observation is available. The strong structure on the missingness pattern is a consequence of the matched nature of the data.

Figure 2 :

 2 Figure 2: Top-left and right panels: generated individual longitudinal data according to a mixed model for n = 100 randomly selected participants prior to dropout (top-left) and after 6 different generated mechanisms of dropout with rate 0.6 (right). Bottom-left: the mean function and the two functional principal components that can be derived from the top-left individual trajectories using a FPCA.

Figure 3 :

 3 Figure3: Relative RMSEs computed on missing test data over 1000 replicates for four missing data scenarios (MCAR, fixed MAR, increasing MAR and increasing MNAR from left to right), two rates of dropout (30% and 60%) and three frequencies of visits (every 3, 2 and 1 year from top to bottom). Each boxplot summarises the RMSEs of a model relative to the reference model (from left to right): LMM quad, LMM cub, LMM spl quant, LMM spl equi, FPCA fve90, FPCA fve99, FPCA k4, JM (only for the MNAR case).

Figure 4 :Figure 5 :

 45 Figure 4: Relative bias over 1000 replicates (left panel), and true plus 10 randomly estimations (right panel) for the mean function and each of the two functional principal components in MCAR and MAR scenarios. The settings are a measure every 2 years for 12 years with a dropout rate of 0 (no dropout), 30% or 60%.

Table 1 :

 1 Index age corresponds to the matching age at which cases were diagnosed with dementia, and controls were still free of dementia. Summary statistics of demographic variables and cognitive scores in the case and control groups from the nested matched case-control (N=348, the 3C Study cohort). Are reported count (frequency) for categorical variables and median (interquartile range) for continuous data and number of measures.

	Variable	Cases	Controls
		(N=174)	(N=174)
	Women	129 (74.1%)	129 (74.1%)
	Index age *	84.8 (81.6 -88.9) 84.8 (81.6 -88.7)
	Secondary school or higher	106 (60.9%)	106 (60.9%)
	At least one APOE-ϵ4 allele	44 (25.3%)	44 (25.3%)
	Follow-up at diagnosis (in years)	9.72 (6.78 -12.2) 9.70 (7.08 -12.4)
	MMS score at index age	23 (21 -25)	28 (25 -29)
	Number of individual measures of MMS	6 (4 -7)	6 (5 -7)
	BVRT score at index age	9 (8 -11)	11 (9 -12)
	Number of individual measures of BVRT	5 (4 -6)	5 (4 -7)
	IST30 score at index age	28 (22 -34)	39 (34 -47)
	Number of individual measures of IST30	5 (4 -7)	6 (4.25 -7)
	TMTA score at index age	0.25 (0.21 -0.35) 0.39 (0.31 -0.48)
	Number of individual measures of TMTA	4 (3 -5)	4 (3 -5)

*

Acknowledgement

This study received financial support from the French government in the framework of the University of Bordeaux's France 2030 program / RRI PHDS. Computer time for this study was provided by the computing facilities of the MCIA (Mésocentre de Calcul Intensif Aquitain).

Additional information

A supplementary material is available, providing details on the simulations scenarios and additional results. Code used for the simulation is shared on the Github repository fpcarobustness.


Jean-Marie Chambard
email: jean-marie.chambard@sanofi.com

Tânia C Gonçalves

Evelyne Benoit

Rachid Boukaiba

Michael Kurtz

Gerhard Hessler

Andrees Bohme

Rémy Béroud

Michel De Waard

Denis Servent

Tânia Goncalves

Rémy Beroud

Michel Partiseti

Venoms in drug discovery: a case study on voltage-gated channels

Keywords: Cyriotoxin-1a, sodium channel, Nav1.7

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.