Catherine Huber-Carol
email: catherine.huber@parisdescartes.fr

Shulamith T Gross
email: shulamith.gross@baruch.cuny.edu

Filia Vonta
email: vonta@math.ntua.gr

Risk analysis: Survival data analysis vs. machine learning. Application to Alzheimer prediction

Keywords: Survival data analysis, Stochastic models, Machine learning, Neural networks, Nonlinear modeling, Alzheimer disease

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.

Introduction

When analyzing the risk of an event E to occur, such as a degradation, a failure, a disease or even death, one may consider how the waiting time Y of onset of such a nocuous event is influenced by intrinsic and environmental factors X := (X 1 , .., X p ):

Y = f (X) (1)
The function f is not deterministic. In classical survival data analysis, a stochastic model for f is chosen among several families of models, fully parametric, nonparametric or semi-parametric [START_REF] Lawless | Statistical models and methods for lifetime data[END_REF][START_REF] Per K Andersen | Statistical models based on counting processes[END_REF]. Then, f is to be estimated, based on the chosen model and from the observation of a training set, i.e. a sample of n observations, (X, Y ) i i = 1, • • • , n, of the pair (X, Y ). Finally, the stochastic model initially chosen has to be tested for fit; see E. J. Lehmann [START_REF] Erich | Testing statistical hypotheses[END_REF], and [START_REF] Huber-Carol | Goodness-of-fit tests and model validity[END_REF]. The classical versions of these models are available in R software. To adapt the analysis to specific situations, researchers have to elaborate extensions of these models and work them out using R, which is both a software and a programming language; see [START_REF] Pons | Estimation in a Cox regression model with a change-point according to a threshold in a covariate[END_REF][START_REF] Terry | Modeling survival data: extending the Cox model[END_REF]. As a counterpart, the machine learning approach of this same problem consists in using an algorithm which has the potential risk factors X as entries and Y as output. Several families of algorithms are available: neural networks, random forests or support vector machines [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF][START_REF] Le Cun | Deep learning[END_REF]. This second approach leads to the so-called "data driven models". In that respect, it seems to be more satisfactory than the subjective choice of a stochastic model that appears in the first approach. However, machine learning is often viewed as a "black box" as the algorithm goes back and forth until convergence is achieved, and it scatters thus the initial potential risk factors in such way that interpretation becomes difficult. Also, every machine learning method, even though it seems to be purely algorithmic, has a probabilistic interpretation. We shall see this feature in particular for neural networks, which are a parametric version of a stochastic model: the projection pursuit regression and discrimination model. Several recurrent problems occur when dealing with duration data: (i) Incomplete data occur very frequently. For some subjects or items, the event does not occur before the end of the study, leading to missing data that are called "censored" data. More precisely, they are called "right censored data" as other phenomenons may occur like left censoring or interval censoring [START_REF] Huber | Interval censored and truncated data: Rate of convergence of NPMLE of the density[END_REF]. Right censoring means that the true value of Y is unknown but known to be "to the right" of the observed duration. Those incomplete data are not to be thrown away. They are taken into account using special devices.

(ii) Specially in medicine, it may happen that the training sample is small. If the sample size is small, it does not help in drawing conclusions. This happens with prospective studies. The simplest example of this case is a clinical trial: only one risk factor, the treatment. It compares a new treatment X = 1 to the usual one X = 0 to increase the life length of patients suffering from a specific disease. For ethical reasons, the sample size n is rather small and the patients are carefully chosen as the characteristics of the two groups should be comparable except for the treatment. This is achieved either by randomization (a double blind study, an experimental procedure in which neither the subjects nor the experimenters know which subjects are in the test and control groups), or pairing on possibly relevant factors [START_REF] Shulamith | Matched pair experiments: Cox and maximum likelihood estimation[END_REF] in order to ensure the comparability of the two sub-samples except for the treatment. For ethical reasons, the size of such samples is necessarily rather small. We should notice that a perfect comparability is achieved in the randomization case only asymptotically, and in the pairing case only if the set of matched factors is correctly chosen. In order to work with an increased sample size, one may use bootstrap techniques [START_REF] Shulamith | Bootstrap methods for truncated and censored data[END_REF], which consist in drawing at random new samples out of the unique sample we have. Several theorems are needed to justify, in each specific case, the use of such duplications of a unique sample; see P. Hall [START_REF] Hall | The bootstrap and Edgeworth expansion[END_REF].

(iii) But when the training set is not carefully collected on purpose like in the preceding example, one has to deal frequently with immense data bases, even in the medical field. Sample size n and/or number of factors p may be very big. It happens when the data are obtained retrospectively, for example from an internet national or international data base. In France, big data sets of patients, with extensive information, are available from the national health insurance. Depending on the goal aimed at, there is a need to extract the useful information, reducing thus the high dimensionality of the data.

The next section, section 2, deals with several types of stochastic survival models. Section 3 gives an illustration of the difficulty of model selection. In section 4, we present some data reduction methods for the high dimensionality of data sets and in section 5 we define neural networks as a parametric version of a stochastic model: projection pursuit regression and discrimination. Finally, in section 6, we give an example of comparison of a stochastic model and a neural network to predict Alzheimer occurrence among patients at Pitié-Salpêtrière Hospital in Paris (France).

Stochastic survival models

The probability distribution of the waiting time Y may be defined by anyone of five different equivalent functions: its survival function S(t) := P (Y ≥ t), its distribution function F (t) := P (Y < t), its density f (t) := F ′ (t) := -S ′ (t) (whenever it is assumed to exist, which is mostly the case), its hazard function h(t) := f (t)/S(t) and, finally, its integrated hazard H(t) := t 0 h(s)ds. The relationships between them are

S(t) = 1 -F (t) = exp(-H(t)) h(t) = H ′ (t) ; H(t) = -ln(S(t)) (2)
Most current models are based on the hazard rate h, the probability that the event takes place at time t, knowing that it did not take place before:

   h(t) = f (t) S(t) f (t) = -S ′ (t) (3)
2.1. Parametric survival models (i) Weibull and generalized Weibull [START_REF] Nikulin | A chi-squared test for the generalized power Weibull family for the head-andneck cancer censored data[END_REF] (a) Weibull model One of the most usual simple parametric model is due to Weibull. It has two positive real parameters λ and α, to be estimated from a training set. Weibull model W (λ, α) is defined as As will be seen below, h is assumed to be a function of p real risk factors X j named covariates in medicine and stresses in industrial setting.

h(t|λ, α) = αλ α t α-1 ; (λ, α > 0) ; t ≥ 0 (4) α = 1 h is constant h = λ (no ageing, ≡ exponential model E(λ) as S(t) = exp(-λt)) 0 < α < 1 h is decreasing ∞ ↓ 0 (see Figure 1 left) α > 1 h is increasing 0 ↑ ∞ (see
(b) Generalized Weibull GW (λ, α, γ)

It has a third parameter γ and reads

h(t|λ, α, γ) = λα γ (λt) α-1 {1 + (λt) α } 1/γ-1 ; (λ, α, γ > 0); t ≥ 0 (5)
GW allows multiple hazard shapes; see Figure 2 GW (λ, 1, 1) = E(λ), Exponential Exponentiated Weibull hazard [START_REF] Efron | Logistic regression, survival analysis, and the Kaplan-Meier curve[END_REF][START_REF] Govind | Exponentiated Weibull family for analyzing bathtub failure-rate data[END_REF] reads

GW (λ, α, 1) = W (λ, α), Weibull If α > 1, α > γ, h : 0 ↑ ∞ If α = 1, γ < 1, h : λ γ ↑ ∞ If 0 < α < 1, α < γ, h : ∞ ↓ 0 If 0 < α < 1, α = γ, h : ∞ ↓ λ t 1 0
h(t|λ, α, γ) = λα{1 -exp[-(λt) α ]} (1-γ)/γ exp[-(λt) α ](λt) α-1 γ{1 -(1 -exp[-(λt) α ])} 1/γ , (λ, α, γ > 0), t ≥ 0 (6)
All moments of EW are finite. Sub-models are EW (λ, α, 1) = W (λ, α) and EW (λ, 1, 1) = E(λ).

For α > 1, α ≥ γ, the hazard h increases from 0 to ∞. For α = 1, γ ≤ 1, the hazard h increases from ( λ γ ) to ∞. For 0 < α < 1, α < γ, the hazard h decreases from ∞ to 0. For 0 < α < 1, α = γ, the hazard h decreases from λ to 0.

(ii) General parametric models

The waiting time Y is defined through a link function g and a random variable Z:

g(Y ) = β T x + σZ (7)
where x are the risk factors, β the parameters to be estimated, g identity or log, and density f of Z may be

             logistic : f (z) = e -z (1 + e -z ) 2 normal : f (z) = 1 √ 2π e -z 2 2 extreme : f (z) = e z-e z (8)
As an example, let us see how Weibull is represented in this setting: if Y is Weibull W (λ, α), g = log and Z ∼ E(1), exponential(1) so that

log(Y ) = β 0 + β T x + σ log(Z) = -log(λ) + 1 α log(Z) (9)
In terms of

W (λ, α) (4): λ = e -(β0+β T x) , α = 1/σ. (iii) Estimation method: Maximum Likelihood (ML) What is known is the observed training set (y i , x i ) i=1,••• ,n
, and a stochastic model

P = L(Y |X),
the probability distribution of Y knowing X, implying unknown parameters to be estimated from the training set.

The unknown parameters are -p parameters:

β = β 1 , • • • , β p quantifying the weight of each risk factor X j , j = 1, • • • , p
-the baseline hazard function h 0 (t) the probability, having "survived" up to time t, to have the event at time t when no risk factors is present. If h 0 is assumed to be a given function of k parameters, like a Weibull depending on k = 2 parameters, we have, in all, p + 2 parameters to estimate.

The values chosen for the parameters are those that maximize the probability of the observed training set, called Maximum Likelihood (ML) estimators. It should be mentioned that one has to be careful when using maximum likelihood looking at the ratio of the parameters to estimate as compared to the size of the training sample which should not exceed some value [START_REF] Huber | Remarques sur le maximum de vraisemblance. Qüestiió: quaderns d'estadística i investigació operativa[END_REF] As the likelihood of one observation, L(α, λ, β|y,

x 1 , x 2 , • • • , x p ) is equal to f (y, α, λ, β, x)
, and all observations are assumed to be independent, the likelihood of the observed training set is equal to the product

L = n i=1 L(α, λ, β|y i , x i ) = n i=1 f (y i , α, λ, β, x i
). Now, find the values of α, λ, β that maximize the likelihood L:

α, λ, β = argmax α,λ,β L(α, λ, β|y, x) (10)
Actually as it is easier to find the derivative of a sum than of a product, in order to maximize an expression, one deals with the log of the likelihood (named log-likelihood) rather than the likelihood itself.

(iv) Now, how to deal with right censored data? Due to the nature of the data, it may happen that Y is not observed when the experiment stops at time C before the event takes place. One can then provide the fact that what is known is that Y > C and replace in the likelihood, the density f at the unknown time Y by the survival S at the known time C. Other types of censoring may occur besides right censoring: left censoring, when the duration is known to be bigger than some observed duration C, interval censoring, when the duration is known to lie between two observed values, C 1 and C 2 . Moreover, what may also happen is truncation: some items or patients may be skipped from the observed sample due to the experimental scheme. In that case, some special procedures have to be applied [START_REF] Shulamith | Regression models for truncated survival data[END_REF][START_REF] Huber | Estimation of density for arbitrarily censored and truncated data[END_REF][START_REF] Huber | Efficient regression estimation under general censoring and truncation[END_REF].

Semi-parametric survival models

Semi-parametric models have a non-parametric part (an "infinite dimensional" unknown parameter: e.g. one or several functions) and a parametric one (a finite dimensional unknown real parameter), both to be estimated through a training set, but focussing on the parametric part, for the sake of interpretation.

(i) Cox model [START_REF] David R Cox | Regression models and life-tables[END_REF][START_REF] David R Cox | Analysis of survival data[END_REF][START_REF] Bretagnolle | Effects of omitting covariates in Cox's model for survival data[END_REF] The hazard rate h is assumed to be equal to a baseline hazard h 0 (t) (the nonparametric part of the model), modified by p covariates X = (X 1 , • • • , X p ) whose weights are the parameters β = (β 1 , • • • , β p ) (the parametric part of the model) to be estimated, as well as h 0 which may be any positive function:

h(t|X) = h 0 (t) e β T X (11)
Usually, one is only interested in β and the objective is to get rid of the nuisance infinite dimensional parameter h 0 . This can be done using Cox's partial likelihood (see below). Once the respective weights of the risk factors have been estimated, one can derive an estimation of the function h 0 .

(ii) Estimation method In this semi-parametric case, the baseline h 0 is free to be any function. The likelihood is then replaced by Cox's partial likelihood: at each failure time y i , consider the ratio of the item that fails to all items still present at risk. Those at risk are those who are neither failed nor censored at time y i . This eliminates h 0 , as it has the same value, h 0 (y i ) for everyone, and takes advantage of all items that will be censored later on, in the denominators. Then, once we have estimated the coefficients β of the risk factors, we can estimate h 0 .

(iii) Frailty models Up to now, the observed population is assumed to be homogeneous and all items (or patients) are independent. In order to model a possible inhomogeneity, one can introduce into the model a random effect η, called frailty [START_REF] Vonta | Efficient estimation in a non-proportional hazards model in survival analysis[END_REF][START_REF] Huber | Frailty models for arbitrarily censored and truncated data[END_REF][START_REF] Huber | Semiparametric transformation models for arbitrarily censored and truncated data[END_REF], acting multiplicatively on the hazard rate h(t|x) of an individual with covariate vector x.

h(t|x, η) = η exp(β T x) h 0 (t) (12)
A frailty model may be considered as a Cox model with an unobserved covariate ln(η) whose coefficient is equal to 1 and whose distribution function is known, derived from F η the distribution of η. Flexible distributions are usually chosen for F η , the most common being gamma, but also inverse gaussian and all stable distributions, leading to the model

S(t|x, η) = exp(-η exp(β T x)H 0 (t)) (13)
where H 0 (t) is the baseline cumulative hazard. Thus, as η is not observed, what is available is the survival integrated with respect to η:

S(t|x) = ∞ 0 exp(-u exp(β T x)H 0 (t)) dF η (u) = exp(-G(exp(β T x)H 0 (t))) ( 14
)
where G is -log of the Laplace transform of η distribution function

G(y) = -ln( ∞ 0 exp(-uy)dF η (u) (15)
The simple Cox model is obtained when G is the identity.

Nonparametric survival models

In industrial setting, accelerated failure time models (AFT) [START_REF] Bagdonavicius | Accelerated life models: modeling and statistical analysis[END_REF][START_REF] Bagdonavičius | Goodness-of-fit tests for accelerated life models[END_REF] are more popular than Cox model and covariates are called stresses which are often controlled while the covariates are simply observed: for G a survival function, i.e. a decreasing function from 1 to 0 on R + , and r a positive function on a p-dimensional process X

S(t|X(s), 0 ≤ s ≤ t) = G( t 0 r(X(s))ds) ∀X ∈ X (16)
We can see there a totally nonparametric model, involving two unknown functions, G and r. It has a semi-parametric version when one of the two unknown functions involved, G or r, is replaced by a function depending on a finite number of parameters; and a fully parametric version if both functions are assumed to be known except for one or more real parameters to be estimated. Those replacements should be done with flexible functions able to adopt different shapes. Other examples of nonparametric models may be found in [START_REF] Droesbeke | Société mathématique de France, Association pour la statistique et ses utilisations (France), and Journées d'étude en statistique[END_REF][START_REF] Shulamith | Nonparametric estimation and regression analysis with left-truncated and right-censored data[END_REF]. Let us remark that parametric models may seem to be too coercive as compared to nonparametric ones. One way to remedy this possible defect is to replace the distribution function of a parametric model by a neighborhood defined through a distance on probability distributions like Hellinger or Prokhorov distance. This leads to a fourth type of model, called robust models [START_REF] Peter | Robust Statistics[END_REF][START_REF] Huber | Robust versus nonparametric approaches and survival data analysis[END_REF].

Latent variable model: First Hitting Time

Health L(t) of a patient, or operational state L(t) of a technological material is a latent (not observed) variable that decreases to 0 due to three types of risk factors. The event occurs when this process reaches a boundary, usually 0 [START_REF] Ting | Threshold regression for survival analysis: modeling event times by a stochastic process reaching a boundary[END_REF][START_REF] Lee | Threshold regression for survival data with time-varying covariates[END_REF]. In order to make clear the structure of this model, let us give a motivating medical example: One has to estimate the expected years of life free of lung cancer lost due to occupational exposure to asbestos. Such a study was required based on a french case-control survey [START_REF] Chambaz | Analysis of the effect of occupational exposure to asbestos based on threshold regression modeling of case-control data[END_REF]. The model is as follows:

L(t|h, µ) = ℓ + µt + B(t) (17)
where (i) ℓ > 0, the initial amount of health, is a function of the initial covariates X I : gender, past family disease history, genetic factors,...(Note that in our example, the amount of health is meant with respect to lung cancer occurrence, not with respect to death itself).

(ii) µ < 0, the slope of the process, is a function of X I and also of lifetime covariates X L : smoking and food habits, environment, biological measurements (glycemia, cholesterol, • • • ).

(iii) B(t) is a Brownian motion, the random part of the model.

(iv) Finally, the focus is put on a special risk factor X S (here occupational exposure to asbestos), which accelerates by a function r(t) the time to onset of lung cancer.

The initial amount, ℓ, will depend on x I , the slope, µ, on x I and x L . Usually this dependence is assumed to be linear, which results in a parametric model, with the weights of each risk factor to be determined for ℓ and for µ. Now, given the values of x I and x L of a patient, one can compute the years free of disease that he lost due to his exposure to asbestos by using the acceleration function r, estimated from the training set. In the absence of any asbestos exposure the waiting time to the event (lung cancer occurrence) is

T (ℓ, µ) = inf{t ≥ 0 : L(t|ℓ, µ) ≤ 0}, (18)
T (ℓ, µ), (< ∞ as µ < 0) is inverse Gaussian (Φ(t) = t -∞ 1 √ 2π exp(-u 2 /2) is the standard normal cdf (cumulative distribution function)): F (t|ℓ, µ) = 1 + e -2ℓµ Φ (µt -ℓ)t -1/2 -Φ (µt + ℓ)t -1/2 , (19)
In the presence of asbestos exposure, the time to event is accelerated by a nondecreasing continuous function r on IR + such that r(t) ≥ t ∀t:

T (ℓ, µ, r) = inf{t ≥ 0 : L(r(t)|ℓ, µ) ≤ 0} (20)
This parametric model has a nonparametric version [START_REF] He | First-hitting-time based threshold regression[END_REF].

Illustration of the difficulty of model selection

Model selection is a difficult problem. The following artificial example is an illustration of the difficulty of eliminating irrelevant predictors. It shows that removing potential risk factors which seem to be independent of the outcome may lead to serious errors. Let us consider a very simple diagnosis problem, discrimination between 2 diseases M = M 1 or M = M 2 based on 3 symptoms X j , j = 1, 2, 3, meeting the following distribution:

M = M 1 X 1 = 0 X 1 = 1 ❍ ❍ ❍ ❍ ❍ X 2 X 3 0 1 0 1/4 0 1 0 1/4 ❍ ❍ ❍ ❍ ❍ X 2 X 3 0 1 0 0 1/4 1 1/4 0 M = M 2 ❍ ❍ ❍ ❍ ❍ X 2 X 3 0 1 0 0 1/4 1 1/4 0 ❍ ❍ ❍ ❍ ❍ X 2 X 3 0 1 0 1/4 0 1 0 1/4
(i) All 3 symptoms are present with the same probability 1/2 in M 1 and in M 2 . This implies that all 3 symptoms seem to be independent of the type of the disease, M 1 or M 2 .

(ii) For any pair of symptoms x j , x k , P (M 1 |x j x k ) = P (M 2 |x j x k ) = 1/4, so that every pair (X j , X j ′ ) is uniform on its 4 values in M 1 as well as in M 2 . As a result, none of the 3 pairs (X j , X j ′ ) can discriminate M 1 and M 2 .

(iii) But for any triplet of symptoms (x 1 , x 2 , x 3 ), whenever P (M 1 |x 1 , x 2 , x 3 ) = 1), P (M 2 |x 1 , x 2 , x 3 ) = 0 and the reverse is also true. This implies that (X 1 , X 2 , X 3 ) altogether discriminate perfectly M 1 and M 2 and lead to a perfect diagnosis of M .

This phenomenon can be extended to any k-uple of risk factors, so that one has to take into account the maximum possible size of sets of risk factors allowed by the number of observed items as compared to the number of risk factors.

Reduction methods

Among existing methods (POD, PGD and others [START_REF] Ohayon | Reduced models for fluid-structure interaction problems[END_REF]), let us cite PGD (Proper Generalized Decomposition method) of Chinesta, Ladevèze et al [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF]:

f (x 1 , .., x p ) ≃ q i=1 ( p j=1 f j i (x j )) (21)
and the classical Singular Value Decomposition (SVD) of a matrix A n×p , which is useful when the data is a high dimensional (rectangular) (n × p) matrix A. It replaces A by a product of three matrices, the inner one being diagonal with dimension r, the rank of A, smaller than min(n,p). It can be seen as a Principal Component Analysis of the matrix A n×p , which replaces the explanatory variables X for the outcome Y by a few linear combinations of them, giving thus a probabilistic interpretation of SVD. Finally, I will cite a decomposition of a function f (x 1 , .., x p ) of several discrete variables, based on variance analysis, that is useful for sparse contingency tables. For an overview of existing methods, see [START_REF] Nouy | Low-rank tensor methods for model order reduction[END_REF].

Singular Value Decomposition (SVD)

A : R p → R n rank(A) = r ≤ min(n, p) :

A n×p =        
                       n : n ≪ p, a centered matrix p Find U , V , D such that: A = U DV T or U T AV = D ( 22
)
where

     U n×r orthonormal basis of E r R n column space of A V p×r orthonormal basis of E ′ r R p row space of A D r×r diagonal matrix (d 1 , • • • , d r ) ; r = rank(A) (23)
Let us remark that U and V have left inverses:

V T r×p V p×r = U T r×n U n×r = I r×r (24)
In order to obtain U , V , D verifying Eqs. ( 22), we proceed as follows:

As U T U = I r×r , A T A = V D 2 V T and, as V T V = I r×r , AA T = U D 2 U T . A T A,
                      n : n ≪ p, a centered matrix p
A represents n observations (on n patients or n technological devices) of the p-dimensional variable (covariate or stress) X = (X 1 , • • • , X p ), and it is considered as a cloud of n points in R p or p points in R n . The goal is to find the best representation of the cloud in a space of smaller dimension, while keeping most of the initial information. This suggests to maximize the variance of linear combinations of the vectors X j , which is obtained by projection on the first k eigenvectors of the symmetric semi-definite positive matrix A T A. This matrix is an estimator of Σ, the unknown covariance matrix of X. This can be seen from the following equations: Let W be any linear combination of the

X ′ j s, W = a 1 X 1 + • • • + a p X p It has a representation in terms of the eigenvectors, W = b 1 V 1 +• • •+b r V r , so that V ar(W ) = a T A T Aa = b 2 1 λ 1 +• • •+b 2 r λ r , where r j=1 b 2 j = 1.
We can see from those equations that the "most informative k-dimensional subspace" in R p is the space spanned by (

V 1 , • • • , V k ), the first k eigenvectors of A T A.

Decomposition of a function of several random variables

Let f be a real function of p real random variables. We have the following decomposition lemma based on the classical variance analysis: Every integrable function f = f (x 1 , • • • , x p ) may be decomposed uniquely into a sum:

f = C + p j=1 g j (x j ) + j,k∈{1,••• ,p} 2 , j<k g jk (x j , x k ) + .. + g 1,..,p (x 1 , .., x p ) ( 25
)
where C is a constant and all expectations of g functions on any of their arguments are 0. This result will be essentially applied in the following case: the variables X j are discrete and take a finite number of values. The expectations are taken with respect to the joint probability of X = (X 1 , •, X p ). Constructive decomposition

C = E(f (X 1 , .., X p )) g j (x j ) = E(f (X 1 , .., X p )|X j = x j ) -C g j,k (x j , x k ) = E(f (X 1 , .., X p )|(X j , X k ) = (x j , x k )) -g j (x j ) -g k (x k ) -C etc . . .
This approach is very useful for sparse contingency tables as can be seen from the following diagnosis example, where the data leave many symptom profiles empty. We have p = 9 binary syptoms X j , j = 1, • • • , p, X j = 1 if the symptom is present, which gives p ′ = 2 9 = 512 symptom profiles, which are supposed to discriminate two different forms, M 1 and M 2 , of the same disease M .

Observations on n = 150 patients give counts for 1024 cells, so that many cells are empty in the observed sample, though their probabilities may not be assumed to be equal to 0.

A 2×512 =   n
     q = 2 p ′ =512 log(P (X = x|M 1 )) = C + p j=1 g j (x j ) + j =j ′ , j<j ′ g j,j ′ (x j , x j ′ ) + j =j ′ =k, j<j ′ <k g j,j ′ ,k (x j , x j ′ , x k ) + • • • + g 1,2,••• , p (x 1 , x 2 , • • • , x p )
We have two such developments, one for disease M 1 , and one for disease M 2 , with sums of functions of an increasing number, k, of variables. If we can assume that the symptoms are independent, which is very simple but often unrealistic, we can stop the preceding developments at k = 1. If we stop at k = 2, it means that we assume a possible dependence of order 2, but with no influence of any third symptom on the link between any two given symptoms. Now, stopping at k, named assumption H k , i.e. cutting off all functions of more than k arguments, means that we allow possible order k interactions, but no interaction greater than k. Under assumption H k , the marginals of order k of the data are sufficient statistics [START_REF] Huber | Estimation dans les tableaux de contingence a un grand nombre d'entrées[END_REF].

Neural Networks

The most usual Neural Network, the single hidden layer back-propagation network (or single layer perceptron), is a particular case of a stochastic model, the Projection Pursuit Regression and Discrimination (PPRD) model.

Projection Pursuit Regression and Discrimination model (PPRD)

This semi-parametric model is fit for regression as well as discrimination problems.

(i) Regression

The target Y ∈ R is the response variable to X = (X 1 , • • • , X p ) ∈ R p . The PPR Y of Y is defined as:

Y = f (X) := M m=1 g m ( ω T m X) := M m=1 g m (V m ) (26
x i , y i ), i = 1, • • • , n.
This is an additive model, but not with respect to the initial variables X but with respect to appropriate linear combinations of them: V m = ω T m X. A remarkable property of this model is that if we choose M big enough, any continuous function may be approximated arbitrarily well. This is what made the success of neural networks which are a parametric version of this stochastic model. A simple example of the nonlinearity of the model may be given in the simplest case where p = M = 2,

ω 1 = (1/ √ 2, 1/ √ 2), ω 2 = (1/ √ 2, -1/ √ 2), g 1 (t) = t 2 /4, g 2 (t) = -t 2 /4, so that the resulting value of f is f (X) = X 1 X 2 .
However, a drawback of this model is the difficulty of interpretation of the results in terms of the initial risk factors as each feature X j is scattered into every linear combination of X. The quadratic measurement error is:

R(θ) := n i=1 [y i - M m=1 g m (ω T m x i )] 2 (27)
where θ is the set of parameters of the problem (θ := (ω m , g m )).

(ii) Discrimination into K categories Up to now we considered the case of a regression problem, Y being a random variable in R. If the problem is a discrimination one, the target Y = (Y 1 , • • • , Y K ) is one of K categories, each Y k being coded as a (0,1) variable. In that case, two different error measurements are considered:

R 2 (θ) := K k=1 n i=1 (y ik -f k (x i )) 2 quadratic error R KL (θ) := -n i=1 K k=1 y ik log(f k (x i
)) crossed entropy The index KL for the crossed entropy makes reference to Kullback Leibler distance which is equal up to an additive constant to crossed entropy as defined up there. We recall that the Kullback-Leibler distance of two probabilities P and Q is defined as

KL(P, Q) = log( dP dQ ) dP 5.1.1.
Neural Network as a special case of PPRD Let the framework be a discrimination problem: the target Y is a category, each Y k being coded as a (0,1) variable. Y k is modeled as a function g k of a linear combination of variables obtained by a linear combination of activated M linear combinations of the inputs. We can see that linearity comes in twice in this definition, with p × M coefficients α and p × K coefficients β. To those parameters are added the activation function σ and the K functions g k :

V m := ω T m X := α 0 + α T m X m = 1, 2, • • • , M Z m = σ(V m ) σ is the activation function T k = β 0k + β T k Z k = 1, 2, • • • , K f k (X) = g k (T ), k = 1, 2, • • • , K
where g k (t) = Average correct predictions due to GLM and NN are p d for dementia cases, p nd for people without dementia, p g for global and CI 95% are the respective confidence intervals.

Hereafter is an illustration, see Figure 4, of the neural network we used, with 2 layers, the first one with 3 neurons and the second one with 2 neurons. In order to have a more readable picture, we choose to show as entries only the four relevant predictive factors: age, incapacity, depression and the gene APOE4, a reduced number (only 4) as compared to the initial 17 entries. Some comments (i) It was very surprising that the education level did not play a role in the occurrence of the disease. (ii) The fact that the counts are very much unbalanced (142 dementia cases versus more than 4000 without dementia) creates problems for the prediction: the confidence intervals are large. (iii) To overcome this problem, one can duplicate the smaller category [START_REF] Le | Personal communication[END_REF], see (iv) The prediction abilities of the two methods are very comparable on this example, with moderate dimension of the data set.

(v) Increasing the number of neurons from 2 to 3 in a unique hidden layer does not improve the prediction ability. The only consequence is increasing the time to get the result. Same comment when using two layers with respective number of neurons 3 and 2

As a perspective, one could use world wide data (Big Data) related to Alzheimer disease and the observed risk factors on the patients under survey. Then this would become a regression problem, estimating the time to onset of the disease as a function of the numerous available risk factors.

Conclusions and perspectives

(i) Analysis of risk is a crucial issue nowadays. The same mathematical methods to evaluate the impact of risk factors on the waiting time for an event to occur may be applied to technological aging systems [START_REF] Rykov | Reliability of engineering systems and technological Risk[END_REF], to public health problems [START_REF] Vonta | Statistical models and methods for biomedical and technical systems[END_REF], insurance, management [START_REF] Harlamov | Stochastic risk analysis and management[END_REF], ecology and other fields.

(ii) The link between survival analysis and reliability, ignored for a long time, so that the terminology is different in both fields to name the same concepts, is now fully acknowledged. Only recently, their common points, more numerous than their differences, were recognized. The favorite model though for survival data is the Cox model and its extensions, while accelerated models are preferred in an industrial environment [START_REF] Bagdonavicius | Accelerated life models: modeling and statistical analysis[END_REF][START_REF] Bagdonavičius | Goodness-of-fit tests for accelerated life models[END_REF].

(iii) The increasing power of computers facilitates the development of non parametric and semi-parametric stochastic models, more greedy in time computation than the parametric ones, and also of complex algorithms in Artificial Intelligence to deal with Big Data. Many new packages appear every month in R and Python which are available on the web.

(iv) Statistical Learning and Big Data (SLBD):

Research statisticians turn now to Machine Learning, like Hastie, Tibshirani et al [START_REF] Hastie | The elements of statistical learning: data mining, inference and prediction[END_REF] and the reverse is also true. Machine learning researchers like Murphy [START_REF] Murphy | Machine learning: a probabilistic perspective[END_REF] investigate the probability background of their manipulation of Big Data through algorithms. This forms finally what could be called Statistical Learning and Big Data.

(v) The increasing knowledge of the genome [START_REF] Ting | Analysis of microarray gene expression data[END_REF] provides high dimensional data [START_REF] Ting | Analysis of microarray gene expression data[END_REF]. The FHT (First Hitting Time) model, in use in this context, gives a special status to some risk factor.

(vi) Specific developments in medicine: "Individualized medicine", which seems to be in total contradiction with classical statistical analysis is now seriously taken into account. Including all the data relative to a patient, such as medical doctors report (text data) leads to the development of new graphical models.

Figure 1

 1 right).

Figure 1 .

 1 Figure 1. Weibull hazards for diverse values of α (shape), λ (scale)

Figure 2 .

 2 Figure 2. Generalized Weibull hazards for diverse values of the parameters

4 . 2 .x

 42 the Gram matrix associated to A, and AA T are semi-definite positive matrices, so that the columns of V are the real eigenvectors of A T A and the columns of U are the real eigenvectors of AA T and d 2 1 , • • • , d 2 r are the common positive eigenvalues of A T A and AA T . As a consequence, the singular value decomposition of A is U DV T . Statistical approach to SVD: Principal Component Analysis (PCA) X 1 X 2 ....... ...... X p 11 x 12 ............... x 1p x 21 x 22 ............... x 2p ... ... ............... ... x n1 x n2 ............... x np 

  ) in which ω m , m = 1, • • • , M are M unitary p-dimensional vectors and functions g m : R → R. Function g m (V m ) is called a ridge function. The estimators are based on the observed training sample: (

e

  T i . This choice ensures that all Y k := f k (X) are positive and add to 1. There are several possible choices for the activation function σ, see Figure3. All of them are smoothed versions of the step function s(u) = 1 {u ≥ 0} (up to an additive constant). The nonlinearity of the model is due to the activation function. If σ is the identity, the model becomes linear.σ(u)= 1 1 + e -u the sigmoïd, the most usual one σ(u) = e u -e -u e u + e -u hyperbolic tangent (th(u))σ(a, u) = a(e u -1) for u < 0 u for u ≥ 0 Exponential Linear Unit (ELU) σ(a, u) = au for u < 0 u for u ≥ 0 Rectified Linear Unit (ReLU) σ(a, b, u) = ba(e u -1) for u < 0 u for u ≥ 0 Scaled Exponential Linear Unit (SELU)

Figure 4 .

 4 Figure 4. 2 layers, 3 and 2 neurons

  11 n 12 ............... n 1p ′ n 21 n 22 ............... n 2p ′

Table 1 :

 1 Method p d p nd pg CI 95% (p d ) CI 95% (p nd ) Comparison of prediction abilities of GLM and NN.

	GLM 0.72 0.73 0.73 0.55 0.85 0.70 0.76
	NN(2) 0.68 0.73 0.73 0.50 0.85 0.65 0.77

Table 2 :

 2 Method p d p nd pg CI 95% (p d ) CI 95% (p nd ))

	GLM 0.73 0.73 0.73 0.71 0.76 0.71 0.75
	NN(2) 0.75 0.72 0.73 0.73 0.78 0.70 0.75

Table 2 :

 2 correct predictions due to GLM and NN for dementia cases (p d ), for people without dementia (p nd ), global p g , and 95% confidence intervals after duplication.After duplication, the widths of the 95% confidence intervals are reduced [0.71 0.76] instead of [0.50 0.85] for the future Alzheimer detection and [0.70 0.75] instead of [0.65 0.77] for the future non Alzheimer.

Acknowledgements

The authors dedicate this work to professor Mikhail Nikulin, from Saint Petersburg, Russia, who shared our research activities during many years.

Application to Alzheimer prediction

We analyze the comparison of prediction abilities of a stochastic model, the logistic regression (GLM), and a neural network (NN) approach, for the prediction for patient to develop Alzheimer in the next 4 years , based on the observation of 17 risk factors. The data set is issued from Pitié-Salpêtrière Hospital in Paris (France). It has n = 5003 patients each with p = 17 covariates, including 3 genetic factors, familial disease anteriority and personal factors (age, sex, education level, • • • ). 467 patients were excluded for missing values, 142 developed the disease and 4214 were free of disease, and thus called controls [START_REF] Huber-Carol | Within the sample comparison of prediction performance of models and submodels: Application to Alzheimer's disease[END_REF]. Let us notice that this is a very simple problem, which deals only with a discrimination problem between two outcomes (Alzheimer or not) and a data set which is not big at all as it involves only a (5003 × 17) matrix. Note also the very unbalanced counts for diseased (142) and controls (4214) which makes the discrimination uneasy for both methods. The logistic model reads

We proceed as follows (i) Split at random 3/4 of the data set to be the training set. This is a rather huge training set due to the paucity of patients who developed the disease.

(ii) The remnant 1/4 will be the test set on which to predict who will be Alzheimer.

(iii) Use separately logistic model (GLM) and neural network (NN) on the training set to estimate the probabilities R GLM (X) and R N N (X) to develop Alzheimer based on the risk factors X.

(iv) Predict, on the test set, who will be Alzheimer based on the estimations done with both methods. The result is four counts for each method: true positive, false positive, true negative, false negative. Or, equivalently, the probabilities of correct classifications of Alzheimer (p d ) and non Alzheimer (p nd ) (v) Repeat this process N times, for both methods, to obtain confidence intervals for the probabilities of correct prediction of Alzheimer, see Table 1:


Azaddien Zarrouk
email: azaddien.zarrouk@insa-cvl.fr

Karim Belharet
email: karim.belharet@yncrea.fr

Omar Tahri
email: omar.tahri@insa-cvl.fr

Vision-based magnetic actuator positioning for wireless control of microrobots

This work is concerned with targeted drug delivery inside the human body using magnetic microrobots. It proposes a vision-based magnetic platform for guiding microrobots in both open-loop/closed-loop schemes. The open-loop scheme can be used for example in the case of the inner ear, where the microrobots cannot be localized in real time. On the other hand, for more accuracy, closed-loop scheme can be used for organs as the human eye since microrobots can be localized using a vision sensor. For both schemes, the platform is designed to compensate for human body movements. It is composed of a new magnetic actuator mounted on a robot end-effector and a hybrid vision system. The latter consists of a camera and two microscopes, while the newly proposed magnetic actuator is built using four permanent magnets. The proposed actuator has been designed to create a local maximum of the magnetic field magnitude in a planar workspace. This results in a convergence point for magnetic microrobots that are in its influence zone, making possible open-loop control with a satisfactory accuracy.

The procedures for calibrating each component of the proposed platform are described and validated. Finally, several experiments have been carried out to validate the modeling part and to show the feasibility of the concept. The obtained experimental results

show that using such platform, the microrobots guiding can be achieved in open-loop under reasonable perturbations and in closed-loop with an accuracy of 200µm.

INTRODUCTION

Manipulation of objects using an external magnetic field inside human body is a strong challenge in the drug delivery applications. Therapeutic and diagnostic possibilities offered by such technology are growing with the clinical need to make procedures less invasive for the patients. There are many previous works discussing manipulations of magnetic particles using an external magnetic field for medical use. Ophthalmologists were the first to explore manipulation of magnetic objects in 1842. The first step in intravascular magnetic guiding was accomplished in the early 1950 [START_REF] Tillander | Magnetic guidance of a catheter with articulated steel tip[END_REF]. The use of magnetic catheters to remove a foreign body in the stomach appeared in the 1940s [START_REF] Equen | Magnetic removal of foreign bodies from the esophagus, stomach, and duodenum[END_REF] [START_REF] Gillies | Magnetic manipulation instrumentation for medical physics research[END_REF]. More recently, magnetic propulsion of capsule endoscopy has been also under investigation [START_REF] Chiba | Moving of a magnetic actuator for a capsule endoscope in the intestine of a pig[END_REF][START_REF] Morita | In vivo trial of a driving system for a self-propelling capsule endoscope using a magnetic field (with video)[END_REF][START_REF] Kósa | Swimming capsule endoscope using static and rf magnetic field of mri for propulsion[END_REF]. Magnetic guidance was also used in neurosurgery, especially for steering of magnetic catheters into cerebral vessels [START_REF] Montgomery | Superconducting magnet system for intravascular navigation[END_REF].

For all applications mentioned above a magnetic field is used to produce force or torques on a magnetic object located inside organ of the human body. To make a magnetic object reach its desired location, the external magnetic field has to be shaped by an accurate control currents or/and positioning of electromagnets or permanent magnets. The magnetic field of an electromagnet is produced by an electric current unlike the permanent magnets that need no power and no cooling device. This makes it possible to design a portable and light magnetic actuator based permanent magnets, which can mounted on the end-effector of robotic manipulator. However, in practice it is much more challenging to design a permanent-magnet based actuator with a well defined and accurate shape of the magnetic field. As related work, Mahoney et al [START_REF] Mahoney | 5-dof manipulation of an untethered magnetic device in fluid using a single permanent magnet[END_REF] developed a control of a magnetic capsule endoscope, navigating in fluid by using a single permanent magnet coupled to a robotic manipulator. In [START_REF] Mahoney | Generating rotating magnetic fields with a single permanent magnet for propulsion of untethered magnetic devices in a lumen[END_REF], the same authors demonstrated the ability to control an untethered magnetic microrobot using a single rotating permanent magnet. In [START_REF] Simi | Design, fabrication, and testing of a capsule with hybrid locomotion for gastrointestinal tract exploration[END_REF][START_REF] Valdastri | A magnetic internal mechanism for precise orientation of the camera in wireless endoluminal applications[END_REF][START_REF] Valdastri | A novel magnetic actuation system for miniature swimming robots[END_REF], a large magnet was displaced and oriented either manually or using a robotic platform to guide different versions of capsule endoscopy. A comparison between robotic versus manual magnetic steering is proposed by the same authors [START_REF] Ciuti | Robotic versus manual control in magnetic steering of an endoscopic capsule[END_REF]. In the clinical field, the wireless Stereotaxis Niobe magnetic navigation system uses two permanent magnets to guide magnetically tips catheters in the cardiovascular system [START_REF] Carpi | Stereotaxis niobe R magnetic navigation system for endocardial catheter ablation and gastrointestinal capsule endoscopy[END_REF]. The majority of magnetic actuators existing to date are based on closed-loop control for precise motion in the presence of perturbations generated by thermal noise or drifting due to boundary effects. It means that a localization system to recover the position of the microrobot is required. The localization system is very often a vision system, which is not always possible in a closed environment like a human body.In this paper, we propose a magnetic platform, capable of controlling a magnetic microrobot in both an open loop or closed-loop schemes. The first contribution of this work is a four permanent magnet based actuator capable of wireless position control of an untethered microrobot in 2D workspace even in the presence of reasonable perturbations. Its particularity is the creation of local maximum of magnetic field magnitude at a distance that acts as a point of convergence, and allows open-loop control of the magnetic microrobots. The second contribution is the calibration method of a hybrid vision system composed of a camera and two microscopes. This system creates a link between micro scale environment where the microrobots are guided and the macro scale where the magnetic actuator is controlled. There are some techniques in the literature dealing with the calibration of non-overlapping cameras. For example, in [START_REF] Lébraly | Fast calibration of embedded non-overlapping cameras[END_REF], the authors estimate the relative pose between two rigidly linked cameras using estimated motion between two different positions and more. In [START_REF] Kumar | Simple calibration of nonoverlapping cameras with a mirror[END_REF], a mirror is used to overcome the non-overlapping fields of view. The solutions proposed in [START_REF] Lébraly | Fast calibration of embedded non-overlapping cameras[END_REF][START_REF] Kumar | Simple calibration of nonoverlapping cameras with a mirror[END_REF] cannot be applied to our system. First, because our platform is fixed, then the microscopes and the camera are observing objects with different scale. Therefore, using the motion of the vision system or a mirror to observe a calibration pattern by the camera and the microscopes at the same time is not feasible.

The remaining of this paper is organized as follows: the next Section explains the principle and modeling of the proposed actuator. In Section 3, a method for calibrating the hybrid vision system composed of a camera and two microscopes is provided.

Then, Section 4 shows the experimental results using the proposed platform. The paper is concluded in Section 5.

L a s e r

Four permanent magnet based actuator

In the following of this section, it will be shown that with a combination of four permanent magnets a local maximum of the magnetic field can be created at some distance. The obtained local maximum will define a local convergence point inside its influence zone thanks to the produced attractive force by magnetic field gradient.

The magnetic field vectors and their magnitude generated by a system of two permanent magnets in xz-plane at y = 0 are shown in Fig. 2 (a). If the magnetization of the magnets is inverted, the magnetic field vectors change the direction while keeping the same magnitude. The resulting magnetic field combining two permanent magnet systems is the sum of the magnetic fields generated by each the two systems. From this, an idea to maximize locally the magnetic field strength becomes possible if we superimpose the two systems together (Fig. 2 (c) and (d)), in such way the magnets of the second system do not touch physically the magnets of the first system. At the same time, keep their magnetic field strength close to magnetic field strength of the first system to avoid the difference in field strength between the two systems.The magnetic field generated by one rectangular permanent magnet (i) at point x = [x, y, z] T as given by [START_REF] Ravaud | Magnetic field produced by a parallelepipedic magnet of various and uniform polarization[END_REF] is a function of the following parameters:

B i (x) = f (x, p, ω, M mag ) (1)
where p = [p x , p y , p z ] T and ω = [θ , φ , ψ] T are the position and the orientation vectors of the magnet, and M mag is the magnetization of the magnet.

The additive property of magnetic field gives the following equation:

B total = B sys1 + B sys2 (2)
where B total is the resulting magnetic field of the two systems,

B sys1 = B mag1 + B mag2
is the magnetic field of the two magnets of system 1 and B sys2 = B mag3 + B mag4 is the magnetic field of the two magnets of the system 2. As the force is directed from region of low to high magnetic field strength, we will focus on the magnetic field strength to have an idea about the resulting magnetic force directions. From equation ( 2), we can write:

B total 2 = B sys1 2 + B sys2 2 + 2 B sys1 • B sys2 cos(Θ)
where Θ is the angle between the magnetic field of the two systems at point (x, y, z)

in space projected in xz-plane at y=0 (see Fig. 3). From the equation (3), we see clearly that the magnetic field strength at point (x, y, z) in space is maximum when the angle Θ = 0, and minimum when Θ = π. In our case the magnitude of magnetic field generated by the two systems is of the same order, so the angle varies clearly the magnitude of the resulting magnetic field. We obtain the local maximum by fixing the magnets of the system 1, then varying the position and the orientation of the magnets of system 2 while keeping symmetry with respect to y-axis.

The generation of local maximum is shown in Fig. 3, where the magnetic field strength at each point is given by vector magnitude. It is clear that moving away from point [0, 0, 5.93] T × cm in this case, the angle Θ between the magnetic field vectors of system 1 and system 2 varies from 0 until reaching the maximum value of π. Consequently, the magnitude of the resulting vectors decreases in turn. We are interested in the area above z = 5.93 × cm, because we know that when moving away from the magnets, the magnetic field decreases naturally whatever the configuration of the four magnets. Below z = 5.5 × cm, the angle Θ is close to π which minimizes locally the magnetic fields. Fig. 4 shows the magnetic force direction and strength in the area around the convergence point. The trajectory of microrobots placed in different initial positions is shown in this area to demonstrate that without any control efforts, the microrobots converged to a known position, which is the actuator convergence point. visual servoing can be used in such case for more accurate control.

x * 1 0 -2 [ m] x * 1 0 -2 [ m] x * 1 0 -2 [ m] d z * 1 0 -2 [ m] z * 1 0 -2 [ m] z * 1 0 -2 [ m]

Hybrid vision system

Calibration challenges and notations

The vision system is composed of a camera and stereo system of two microscopes.

It has been proposed to meet several needs in the proposed platform, including the need to validate the model of microrobots inside viscous environment, to validate the actuator magnetic model in order to obtain the best performance in the open-loop scheme and finally to control microrobots in closed-loop when they can be localized in real time. In order to use the proposed system, it is necessary to calibrate its different components. Calibrating the hybrid vision system (Fig. 5) consists in estimating the intrinsic parameters of each sensor and the relative poses of the two microscopes with respect to the camera. This has to overpass two challenges. The first is the scale difference between the camera and the microscopes. The second is the non-overlapping fields of the sensors. In the following, our proposed calibration method is explained.

The following notations are adopted for the different parameters to be estimated:

• K cam , K mic1 and K mic2 are respectively the matrices of the intrinsic parameters of the camera and those of the two microscopes.

• cam T mic1 , cam T mic2 and mic1 T mic2 are respectively the homogeneous transformations between microscope 1 and camera, microscope 2 and camera, and microscope 2 and microscope 1.

Recall that the matrix of the intrinsic parameters K is defined by:

K =      α 0 u 0 0 β v 0 0 0 1      (3)
where α and β are the scale factors of the image, u 0 and v 0 are the coordinates of the optical center of the camera in the image frame.

The homogeneous transformation b T a is defined from rotation matrix b R a and translation vector b t a as follows:

b T a =   b R a b t a 0 1   (4)

Calibration method

The main problem in calibrating the hybrid vision system is to estimate the homogeneous transformation cam T mic1 and cam T mic2 , between the camera and the microscopes. To calibrate such vision system, firstly, the camera and the microscopes are calibrated separately in order to obtain their intrinsic parameters. Secondly, we created a rigid link between the microscopes and camera field of view using two patterns rigidly linked as shown in Fig. 6.

The calibration process is as follows: firstly, the calibration pattern is positioned such that the microscopes are able to observe the small grid and the big grid is visible for the camera. For each position of the calibration pattern, images are acquired by each sensor. Then, using a conventional calibration procedure for each acquired image, the corresponding pose of the big grid with respect to the camera cam T M can be recovered.

The poses mic1 T m and mic2 T m of the small grid with respect to the microscopes can be estimated in the same way.

The last step uses the fact that the homogeneous transformation M T m between the two planar patterns is fixed and known to obtain the unknown relative pose between the camera and the two microscopes.

For the relative pose between the camera and microscope 1, decomposing cam T mic1

gives:

cam T mic1 = cam T M M T m m T mic1 (5)
Equation ( 5) can be rewritten as follows:

cam T mic1 = cam T M M T m mic1 T m -1 (6)
Likewise, the homogeneous transformation between camera and microscope 2 coordinate systems can be written as follows:

cam T mic2 = cam T M M T m mic2 T m -1 (7)
Using equations ( 6) and [START_REF] Montgomery | Superconducting magnet system for intravascular navigation[END_REF], the homogeneous transformation between the camera and each microscope can be obtained from a single position of the calibration pattern.

However, for robustness purpose, more than one pattern positions should be used.

Calibration refinement using bundle adjustment

In this part, a procedure to refine the hybrid system parameters using bundle is proposed. For that, it is necessary to define a cost function taking into account reprojection errors in the camera image plane and microscopes image plane as well. Actually, reprojection errors are defined from the difference between extracted points in the image plane and those computed based on the vision sensor model and pattern/sensor pose. In our case, it is necessary to define a cost function considering the microscopes/camera pose in the minimization process, while using a minimal number of parameters. To respect these conditions, we defined three cost functions as follows:

• For the camera:

e cam = n ∑ i=1 m c ∑ j=1 x c i j -f 1 (K cam , cam T M i , X M j ) 2 (8)
• For the microscope 1:

e mic1 = n ∑ i=1 m m ∑ j=1 x m1 i j -f 2 (K mic1 , mic1 T cam , cam T M i , X m j ) 2 (9)
• For the microscope 2:

e mic2 = n ∑ i=1 m m ∑ j=1 x m2 i j -f 3 (K mic2 , mic2 T cam , cam T M i , X m j ) 2 (10)
where:

• n is the number of the acquired images, m c is the number of points extracted from the image acquired by the camera, while m m is the number of those extracted from the microscope image.

• X M j = [x c j , y c j , 0] T is coordinate vector of point j in the big grid coordinate system.

• X m j = [x j , y j , 0] T is coordinate vector of point j in the big grid coordinate system.

• x c i j = [u c i j , v c i j ] T and f 1 are respectively the coordinate of measured point j of image i and the projection of point X wc j in image i acquired by the camera

• x m1 i j = [u m1 i j , v m1 i j ] T , x m2 i j = [u m2 i j , v m2 i j ] T , f
Finally, Levenberg-Marquardt algorithm [START_REF] Moré | The levenberg-marquardt algorithm: implementation and theory[END_REF] is used to solve the nonlinear minimization problem of the cost function e.

Vision based actuator positioning and microrobot steering

This section explains the full process for magnetic microrobot steering by accurately positioning the four-magnet actuator with respect to the micro-workspace using the hybrid vision system. This includes the microrobot trajectory planning in the micro-workspace and the necessary trajectory to be followed by the robot end-effector to position adequately the four-magnet actuator.

Vision-based actuator positioning and microrobot steering strategy

The full process is shown on Fig. 9. First, the trajectory to be followed by the microrobot in the micro-scale workspace is defined and discretized into m intermediate

goal positions m x d = [x d , y d , 0] T .
The coordinates of the points m x d are then expressed in the end-effector frame and noted e f f x d by:

e f f x d = e f f T m m x d (12)
The homogeneous transformation e f f T m can be decomposed as follows:

e f f T m = e f f T target target T cam cam T mic mic T m (13)
where as shown in Fig. 7:

• mic T m : the relative pose of the micro-workspace with respect to the microscopes frame. This later is defined with the aid of markers seen detected with microscopes.

• e f f T target : the relative pose of the target frame with respect to the robot endeffector frame. This rigid pose obtained after robot-camera calibration which a classical problem in robotics field that have been sufficiently studied [START_REF] Lenz | Calibrating a cartesian robot with eye-on-hand configuration independent of eye-to-hand relationship[END_REF][START_REF] Wijesoma | Eye-to-hand coordination for vision-guided robot control applications[END_REF]. • target T cam : the relative pose of the camera with respect to the target obtained thanks to the known structure of the grid.

• cam T mic : the relative pose of the microscopes frame with respect to the camera, obtained thanks to the proposed calibration method of the hybrid vision system.

Since the convergence point has rigid position with respect to the robot end-effector, the microrobot trajectory is expressed at the end-effector because the desired trajectory will be realized by the convergence point, thus by the robot end-effector (see Fig. 8).

The motion to be applied to the robot end-effector is defined by the translation between the convergence point expressed in the end-effector frame and a desired point expressed in the same frame. The convergence point position with respect to the effector may be obtained using the actuator 3D model. In the experimental platform, errors may come from the 3D printing of the actuator and when mounting the actuator on the effector. A calibration step has been performed to avoid these errors thanks to the calibrated hybrid vision system. The calibration step was performed as follows: the actuator mounted on robot manipulator has been driven to position the convergence point in the workspace under the microscopes where a microrobot is placed. As the microrobot converges to this point, its position can be recovered thanks to the microscopic stereo system. Let mic x cp be the position of the convergence point under the microscopic stereo system.

This position is transformed into the robot end-effector using the camera. In fact, the camera creates a link between the frame of the robot and that of the microscopes thanks to the already calibrated hybrid vision system and the robot-camera calibration step. The last calibration problem is a classical problem in robotics field that have been sufficiently studied [START_REF] Lenz | Calibrating a cartesian robot with eye-on-hand configuration independent of eye-to-hand relationship[END_REF][START_REF] Wijesoma | Eye-to-hand coordination for vision-guided robot control applications[END_REF]. Thus, the convergence point position with respect to the end-effector e f f x cp is obtained with the following equation:

e f f x cp = e f f T target target T cam cam T mic m x cp (14)
Once the calibration step is done, the translation motion from the actual convergence point position to a desired position is then calculated as follows:

e f f t ie f f 2 = e f f x id -e f f x cp (15)
The calculated translation is given to the robot to move the end-effector from a desired point to another. As the end-effector moves, the other desired points position has to be updated with respect to the new end-effector frame.

The proposed strategy allows two vision-based tracking possibilities. As seen on Fig. 9 camera can be selected to track markers on the workspace. In this case mic T m is calculated as follows:

mic T m = mic T cam cam T M M T m ( 16
)
where cam T M is the pose of the markers in the camera frame, and M T m is the relative pose of the micro-workspace frame with respect to the markers. From equation ( 14), one can see that positioning the actuator with respect to the workspace is ensured in real time with the calibrated hybrid vision system. However, once the planned trajectory is defined, the positioning is achieved in open-loop with no feedback on the microrobot location. Therefore, the hybrid vision system serves only for positioning the actuator in the workspace. This feature is important in medical field to take into account the movements of the patient during the treatment. For example to guide microrobots inside the inner ear (cochlea), the trajectory is planned initially to be released in open-loop. However, once the patient's head moves, the planned trajectory should be corrected with respect to the new position of the head. On the other hand, the platform can track and correct the microrobot position when it is accessible for visual sensors (see Fig. 9), like the case of the eye treatment. In this case, the microrobot position is estimated with the microscopes, the actual position is compared to the desired position and corrected to minimize the error.

Dynamic modeling of magnetic microrobot in viscous environment

The open-loop navigation of a magnetic microrobot is based on its dynamic model. This model allows to predict the time required for a microrobot to move from a point of the planned trajectory to the next one. In our context, the magnetic microrobot is in a viscous environment and subject to several types of forces: magnetic force, the viscous drag force of Stokes, inertia, buoyancy and gravity, thermal kinetics (Brownian motion), double-layer electrical interactions, and van der Waals force. Among those forces, the magnetic and the viscous ones are the most significant, while the others can be neglected [START_REF] Belharet | Control of a magnetic microrobot navigating in microfluidic arterial bifurcations through pulsatile and viscous flow[END_REF][START_REF] Pankhurst | Applications of magnetic nanoparticles in biomedicine[END_REF][START_REF] Belharet | Endovascular navigation of a ferromagnetic microrobot using mri-based predictive control[END_REF].

Using second Newtons law, the microrobot dynamical model can be written as follows:

m p dv p dt = F mag + F d (17)
where m p and v p are the mass and the velocity of the microrobot. F mag and F d are the magnetic force and the drag force respectively. The magnetic force on a spherical microrobot depends on the magnetic field gradient created at its position [START_REF] Okuhata | Delivery of diagnostic agents for magnetic resonance imaging[END_REF]:

F mag = V m M .(∇ B ) ( 18
)
where V m is the volume of the magnetized object, B is the flux density of the applied magnetic field, and M is the magnetization of the microrobot. On the other hand, the magnitude of the fluidic force acting on a moving microrobot is determined by Stokes law:

F d = -6πηrv p ( 19
)
where η is the viscosity of the fluid, r is the radius of the microrobot, and v p is the velocity of the microrobot. Finally, the movement of a spherical microrobot immersed in a fluid and actuated by an external magnetic field, can be described by a system of ordinary differential equations (ODE):

   ẋ = v p ẍ = V m m p M .(∇ B ) -6πηr m p ẋ (20)
Using [START_REF] Wijesoma | Eye-to-hand coordination for vision-guided robot control applications[END_REF], it is possible to estimate the time necessary for a microrobot to reach a desired position.

Experimental validation

In the following of this section, first, the calibration results of the hybrid system are shown and discussed. Then the convergence point is experimentally proven and its coordinates in the end-effector frame estimated using the hybrid vision system. Finally, experimental results of microrobot steering in open-loop and closed-loop are shown.

Calibration of the hybrid vision system

First, the camera and the microscopic stereo system are calibrated separately and then their relative poses defined by cam T mic1 and cam T mic2 are obtained using equations ( 6) and [START_REF] Montgomery | Superconducting magnet system for intravascular navigation[END_REF]. To further refine these results, the bundle adjustment technique is used to reformulate the optimization problem and Levenberg-Marquardt algorithm is used to solve it. The calibration results are evaluated in two different ways. Firstly, using reprojection errors in the image. Then, by testing the ability of the camera to control the position of a Laser beam mounted on a robot end-effector to target points observed by microscopes.

Calibration results before and after refinement

In the following, the extracted points from acquiring images by the two microscopes are compared to those calculated using following reprojection equations:

x m1 = f 1 (K mic1 , mic1 T cam cam T M M T m , X w ) (21)
x m2 = f 2 (K mic2 , mic2 T cam cam T M m T m , X w ) (22)
where x m1 and x m2 are the coordinates of the projected point in the two microscopes image plane and X w is coordinate vector of point in world coordinate system. The transformations cam T mic1 and cam T mic2 have been introduced in the reprojection on images in order to test their accuracy.

For each microscope, the average and the standard deviation of the reprojection errors defined by equations ( 21) and ( 22) have been computed using 20 acquired images with 42 points each. The obtained results are shown in Table 5.1.1. The left column shows the result obtained using equations ( 6) and ( 7) to estimate cam T mic1 and cam T mic2 before refinement. The average reprojection error and the standard deviation in image plane are around 12pixels and 7pixels respectively, which correspond to 634µm and 371µm in the small grid coordinate system. These results may appear huge at first sight but can be explained by the scale difference between the camera and the microscopes.

Actually, small errors in the estimated camera poses correspond to several pixels in the microscope image. Besides, the obtained results can be refined using a global optimization as explained in section 3. The results obtained after refinement are shown in the right column of Table . Compared to the initial results, the average projection error in the image plane is reduced by 92% and 90% for microscope 1 and microscope 2 respectively. In the small grid coordinate system, this error is reduced by 93% for microscope 1 and 91% for microscope 2.

Table 1: Calibration results of the hybrid vision system sensors before and after optimization. hybrid vision sensor, a Laser line beam mounted on the robot end-effector. The process for obtaining the rigid transformations 'Grid 1'/Laser and 'Grid 1'/robot are omitted on purpose and they are assumed to be known. Currently, the goal is to position the Laser beam toward points on the Grid 2 situated on the field of view of the microscopes.

The Laser beam position is controlled using the pose of the grid (refer to Fig. 10) in the camera frame. On the other hand, the 3D coordinates of the points of the Grid 2 (Fig. 10) to be targeted by the Laser are expressed first in microscopes frame and then in the camera frame through the transformation mic T cam . In practice, there exists an infinity of the Laser motion to reach the points under the microscopes. In this experiment, to move the Laser from a position to another, pure translation and pure rotational motions of the robot end-effector has been used. To move from initially targeted point to target another point, the translation vector, rotation angle and axis are calculated using equations ( 23), ( 24) and (25). The obtained results are shown in (Fig. 12), where 9 points have been chosen to be targeted with the two types of end-effector motion presented previously. The Laser is initially targeting point p 0 . The experiments show that the Laser points correctly the chosen 9 points with very small error. Table 2 shows the targeting error both in the image plane (pixel) and in the space (µm) for each of the 9 points. The average targeting error is 322.04µm.

t = p f -p 0 (23)
a = p i × p 1 (24)
θ = acos(p i • p 1 ) (25)

Magnetic microrobot guiding using the vision-based platform

The aim of this experimental part is to highlight the existence of the convergence point and to show the feasibility of microrobot steering using the proposed platform in both open-loop and closed schemes. As it is shown on Fig. 1, the used prototyping platform is composed by a magnetic actuator, camera and two microscopes (hybrid vision system), and robot manipulator to perform the positioning of the actuator. The workspace where the microrobot evolve has a circular shape (Fig. 13) filled with viscous liquid composed of a mixture of 80% water and 20% glycerol with a dynamic viscosity of 17.288×10 -4 Ns/m 2 . The used microrobot in the experiments is approximately of To demonstrate the effect of the local maximum, a small permanent magnet is brought close to the workspace to shift the particle from its stable position to another. Then, the particle is released by moving away the small magnet from the workspace. The trajectories followed by the micro-particles from the different initial positions are shown in Fig. 15. As predicted, the particle returns to its stable position, which is the correspond local maxima of the magnetic field magnitude generated by the actuator. The obtained results also show that particle returns to the local maximum following path inflected the actuator force vectors. Actually, the force map projected on the workspace as shown on Fig. 15 have been obtained based on the actuator model.

The previous experimental results have shown that the actuator is capable of converging a microrobot to a specific point in a plane. In the new set of experiment this ability is tested while the particles have to be moved along a well defined trajectory in an open loop. Two types of trajectories have been considered: trajectory composed of linear segments, and a spiral. For each trajectory, the actuator is moved to its initial position ensuring that the convergence point and the first point of the trajectory coincide. As the convergence area covers all the rectangular workspace, whatever the initial position of the particle, its convergence to the starting point is ensured. This capacity is important and could be critical when guiding microrobots inside the human body.

Actually, based on this property, moving the convergence point along the predefined trajectory makes the particle follow the same path. Figures (Fig. the plots, it can be seen that initially, the particle follows the desired trajectory with error less than 100µm for the rectangular trajectory and less than 400µm for the spiral trajectory. At the moment of disturbance at 25 th point for the rectangular trajectory and 90 th point for the spiral trajectory, we observe the peaks that represent the maximal error before adjusting the actuator position with respect to the new workspace position.

After adjusting the actuator position, the error is minimized and the particle has come back to follow the desired trajectory.

The second experiment of this part aims to explore the possibility of using visual feedback on the position of the microrobot when it is accessible. Such scenario can be considered for the case of the eye where its interior can be seen with a vision system.

On the contrary of the open-loop scheme, in the closed-loop control the time needed to follow a trajectory is reduced, because it is no longer needed to calculate the maximum time needed to move the microrobot from one point to another. The obtained results of this experiment are shown in Fig. 18, where a complex trajectory has been considered. In the closed loop, a standard error of 200µm has been tolerated to move to the next point of the trajectory, that's why a small gap is noticed in some points of on the magnetic force generated by the actuator. If the magnetic force is strong enough, the actuator will push a microrobot in ragged or arcuate surface. Otherwise, if the magnetic force is not strong enough, a ragged surface will result in friction that may stop up the robot displacements.

Conclusion

This paper proposed a solution for precise magnetic microrobot trajectory tracking which is a vision-based magnetic platform. The two main components of the platform are permanent magnets based actuator and a hybrid vision system. The magnetic actuator is a four-magnet based actuator capable of generating a convergence point at a distance in the plane, which facilitates the control of microrobots. Its design principle and method to create the local maxima of magnetic field magnitude at a distance are explained. The hybrid vision system is composed by camera and two microscopes used mainly to create a link between a micro scale environment where the microrobots are guided and the macro scale where the magnetic actuator is controlled. This allows to study the dynamics of microrobots in viscous environment and to estimate precisely the convergence point position with respect to end-effector. Microrobots steering strategy is proposed to demonstrate the feasibility of the robotic platform to control the microrobot in both closed environments and others that are accessible for visual feedback.

The platform main components and the steering strategy are validated experimentally.

The results of the proposed calibration method before and after refinement are presented. The existence of convergence point generated by the actuator is demonstrated.

Two experiments using open-loop and closed-loop schemes are carried out to validate the steering strategy using the proposed platform. Future work will be devoted to improving the current actuator by optimizing its structure.

Figure 1 :

 1 Figure 1: Experimental setup shows the different components of the platform.

Figure 2 :

 2 Figure 2: (a) Magnetic field strength and direction of system 1: the heat map represents the magnetic field strength generated by the magnetic system composed of two permanent magnets, positioned as shown in the small scheme at the bottom left. The black vectors represent the normalized magnetic field direction, and the green vector in the small scheme gives the magnetization direction of the magnets. (b) Magnetic field strength and direction of system 2. The difference from the system in (a) is the magnetization direction of the magnets. Magnets in (a) and magnets in (b) have a magnetization in opposite directions. c) and d): Novel four-magnet based actuator superimposing both configurations (a) and (b).

Figure 3 :

 3 Figure 3: Magnetic field vectors in xz-plane of each system (red vector for system 1 and blue vectors for system 2) and resulting magnetic field of the two systems around the created local maximum in xz-plane.

Figure 4 :

 4 Figure 4: The four-magnet system and the area around its convergence point in xz-plane. Black dots are microrobots in different positions, and yellow narrows are the trajectories released by the microrobots because of the effect of the convergence point. The scheme at the bottom-right corner represents the actuator design in the yz-plane.

Figure 5 :

 5 Figure 5: The hybrid vision system with the various unknown homogeneous transformations between sensors coordinate systems.

  2 and f 3 are respectively the coordinate of measured point j of image i and the projections of point X w j in image i acquired by the first and second microscope. Note that differently from f 1 , the homogeneous transformations used in the reprojection f 2 and f 3 were decomposed to include the homogeneous transformation camera/microscopes.The equation (8) minimizes the reprojection error in the image plane of the camera by optimizing its intrinsic and extrinsic parameters. The equations (9) and (10) will minimize the reprojection error in the microscope 1 and microscope 2 respectively by optimizing their intrinsic parameters and relative poses with respect to the camera. The global cost function is defined as the sum of the three previous ones: e = e cam + e mic1 + e mic2

Figure 6 :

 6 Figure 6: Calibration of the hybrid vision system using planar patterns of two different scales: a) transformation decomposition, b) the planar patterns separated by pure translational motions.

Figure 7 :

 7 Figure 7: Different relative poses between the vision based platform components.

Figure 8 :

 8 Figure 8: robot-end effector displacement when moving microrobot from position to another.

Mi c r o

  r o b o t mo v e me n t p l a n n i n g T r a j e c t o r y a d j u s t me n t Wo r k s p a c e Mo v e me n t d e t e c t i o n Ca me r a Ro b o t ma i p u l a t o r A c t u a t o r Cl o s e d -l o o p Cl o s e d -l o o p Op e n -l o o p A c t u a t o r p o s i t i o n i n g Mi c r o r o b o t g u i d i n g Mi c r o s c o p e s

Figure 9 :

 9 Figure 9: Microrobot steering strategy using the hybrid vision system and the four-magnet actuator.

2 5. 1 . 2 .Figure

 212 Figure The proposed platform with Laser mounted on the robot manipulator instead of the magnetic actuator to validate the calibration method.

  As described inFig 11, t is the translation vector to move from p 0 to p 1 ; a and θ are respectively the axis and the angle of rotation to move from p 0 to p 1 ; p f is the projection of point p 0 on the Laser line containing the point p 1 (Fig 11a); p i is the point belong to Laser line containing point p 0 and have the same magnitude of point p 1 (Fig 11b).

Figure 11 :

 11 Figure 11: (a) Calculation method of the pure translation motion. (b) Calculation method of the pure rotation motion.

Figure 12 :

 12 Figure 12: Validation method results: targeting 9 arbitrary points with two types of end-effector motion, pure rotation and pure translation.

Mi c r o r o b o t Ma r k e r sFigure 13 :

 13 Figure 13: The workspace were the microrobot is controlled. The markers placed on the workspace are used by the vision system to detect its movement.

Figure 14 :

 14 Figure 14: Experimental setup: the robot manipulator to position the actuator. The microscope, to only visualize the particle in the workspace, and the workspace is a rectangular shape with dimensions of 2x3cm where the particle is placed.

Figure 15 :P l a n i fi e d t r a j e c t o r y E x p e r i me n t a l t r a j e c t oFigure 16 :

 1516 Figure 15: Convergence point validation results: (a) shows the magnetic particle initial position at a local maximum, (b), (c) and (d) are the particle trajectories when moving back from positions 1, 2 and 3 respectively, to its initial position under the action of magnetic force. The simulation results of the magnetic force direction are printed in the background of the workspace.

16

 16

  (a) and (b)) show the obtained results for the two considered trajectories tracking with the presence of perturbations introduced by bringing a small permanent magnet close to the particle workspace. The obtained plots show the ability of the proposed actuator to achieve magnetic microrobot with an open loop control scheme.5.2.2. Vision-based actuator positioning, perturbation compensation and microrobot steering strategyThe first set of experiments considered trajectory following by a microrobot in static micro workspace.The particles are driven by positioning the actuator using the pose of its rigidly attached pattern with respect to the camera. The following experi-ment considers scenarios of moving micro-workspace. For instance, this can happen in the case of targeted drug delivery a moving internal ear due the patients head movement. Actually, as it is the case for other medical applications, optical trackers fixed on the patient's head can be used to recover this kind of movements. Here, two options are possibles to follow the workspace, either with the camera by attaching markers that are visible by the camera to the workspace. Alternatively, with microscopes by gluing small markers directly onto the workspace. Choosing the second option, two type trajectory types have been released, rectangular and spiral trajectory. Those trajectories have been first defined in coordinates system related to the microscopes, then expressed in the robot end-effector frame taking into account the movements of the micro-workspace. The obtained results of trajectory following are shown in Fig.17

  (a) and Fig.17

  (b), at a moment during the navigation, disturbances are simulated by moving the workspace to assess the trajectory correction using the vision system.Fig.17 (c) and Fig.17

  (d) show the norm of the error between the predefined and the experimental trajectory. ; the error is calculated from start to end point of the trajectory. From

  Mi c r o r o b o t t r a j e c t o r y Mi c r o r o b o t t r a j e c t o r y P e r t u r b a t i o n

Figure 17 :Figure 19 :

 1719 Figure 17: The pictures present the results of the open-loop trajectory tracking and closed-loop workspace position correction, and the error results between the predefined and the experimental trajectory, respectively: (a) and (c) for a rectangular trajectory , (b) and (d) for a spiral trajectory.

Table 2 :

 2 Results of the validation method for the chosen 9 points.

	Points	Targeting error Targeting error
		in image	in space
	Point 0 5.83 pixel	489.3µm
	Point 1 5.38 pixel	430.8µm
	Point 2 7.61 pixel	491.7µm
	Point 3 3.60 pixel	288.8µm
	Point 4 1.41 pixel	188.5µm
	Point 5 3.00 pixel	199.8µm
	Point 6 5.00 pixel	405.8µm
	Point 7 3.16 pixel	237.2µm
	Point 8 1.41 pixel	166.5µm
	spherical form with a radius around 250µm.

Acknowledgments

This work was supported by the CochléRob project funded by Région Centre, HEI campus Centre, and Châteauroux Métropole.

Microrobot navigation in presence of obstacle

The experiment aims to study the effect of obstacles in the microrobot workspace.

The experiment consists in making the microrobot crosses an obstacle formed by "wall of the bristles", while following a linear trajectory back and forth. The goal is to visualize how the microrobt will cross the wall of the bristles. The results of this experiment are shown in Fig. 19 (b),(c),(d),(e) and (f) where the particle initially performs a linear trajectory as desired. Once the particle reaches the obstacles, the magnetic force forces the particle to cross the wall of bristles with a small deviation.

Once the obstacle is crossed, the particle joins its planned trajectory (see Fig. 19 (c)).

The influence on the horizontal navigation due to changes in vertical direction depends