X Bry

C Lavergne

M Tami

EM estimation of a Structural Equation Model

Keywords: EM algorithm, Factor model, Latent Variable, Structural Equation Model

In this work, we propose a new estimation method of a Structural Equation Model. Our method is based on the EM likelihood-maximization algorithm. We show that this method provides estimators, not only of the coefficients of the model, but also of its latent factors. Through a simulation study, we investigate how fast and accurate the method is, and then apply it to real environmental data.

Introduction

Structural Equation Models (SEM) are widely used is as various research domains as psychology, social and behavioral sciences, ecology, chemometrics, etc. A SEM formalizes the interdependence of many Observed numeric Variables (OV) through fewer unobserved ones, referred to as Latent Variables (LV). Every LV is assumed to be underlying a specific set of OVs, which depend on it as well as on extra observed covariates. A SEM is structured through two types of equations, termed measurement equations and structural equations. A measurement equation relates a LV to the corresponding OV's. A structural equation states a hypothesized relationship between LV's. 1 graphs an example of SEM.

Literature widely presents two competing families of methods that deal with SEM's: factor-methods, and component-methods. Among the former family are the classical Factor Analysis, and Jöreskog's SEM estimation technique [START_REF] Jöreskog | A general method for analysis of covariance structures[END_REF] implemented in the LISREL software. These methods use factors as LV's, i.e. variables of which we merely assume to know the distribution (typically standardized normal). They base their estimation on the structure of the covariance matrix of the data according to the model, within a likelihood maximization approach. They estimate all coefficients in the model (linear combination coefficients and variances), but not the values of the factors, which therefore remain unknown. The componentmodel family of methods assumes that every LV is a component, i.e. a linear combination, of its OV's. Note that such a constraint is stronger than the distribution-assumption made on factors. This family includes the classical Principal Component Analysis, Canonical Correlation Analysis and Partial Least Squares (PLS) methods [START_REF] Jöreskog | Recent Developments in Structural Equation Modeling[END_REF], [START_REF] Wangen | A multiblock partial least squares algorithm for investigating complex chemical systems[END_REF], [START_REF] Lohmöller | Latent Variable Path Modeling with Partial Least Squares[END_REF], W. W. [START_REF] Chin | Structural equation modeling analysis with small samples using partial least squares[END_REF], [START_REF] Vinzi | Handbook of Partial Least Squares: Concepts, Methods and Applications[END_REF], but also more recent techniques as Generalized Structured Component Analysis [START_REF] Husson | FactoMineR: An R Package for Multivariate Analysis[END_REF], Generalized Regularized Canonical Correlation Analysis [START_REF] Tenenhaus | Regularized Generalized Canonical Correlation Analysis[END_REF] and THEME Bry et al. (2012), [START_REF] Bry | THEME: Thematic model exploration through multiple co-structure maximization[END_REF]. Factor-methods and PLS-type ones have been compared in several works [START_REF] Jöreskog | Recent Developments in Structural Equation Modeling[END_REF]. The gist is that the latter encounter less convergence problems than the former with small samples. A second advantage is that, since they express every LV as a linear combination of its OV's, and yield the estimated coefficients of that combination, the values of the LV's are estimated, and can also be forecast on new samples, opening the way to easy cross-validation. Such is not the case of factor-methods, which do have yet the theoretical advantage to be based on a proper statistical distributionbased model of data, contrary to PLS-type methods, thus allowing standard statistical tests, which are not possible with the latter. In many areas, it is of essence to be able to estimate the values of LV's on statistical units, since these values allow to efficiently analyze the disparities of units on a reduced number of dimensions. Therefore, we are interested in estimating these values even in the factormodel context. In this work, we adapt the EM algorithm to the SEM estimation problem, in order to get estimates of the factor values. The paper is organized as follows. Section 2 formally introduces the equations of the SEM we deal with. Section 3 applies the EM algorithm to the SEM and derives the estimation formulas. Section 4 first presents a simulation-based study of the performance of the method, with comparison to more classical methods, and then an application to environmental data.

The model

Notations 2.1.1. Data notations

The data consists in blocks of OV's describing the same n units. We consider the following data-matrices and notations:

Y = {y j i }; i ∈ 1, n , j ∈ 1, q Y is the n × q Y matrix
coding the dependent block of OV's y 1 , ..., y q Y , identified with its column-vectors.

X m = {x j,m i }; i ∈ 1, n , j ∈ 1, q m , m ∈ 1, p is the n × q m matrix
coding the m ieth -explanatory block of OV's x 1,m , ..., x qm,m . Value of variable x j,m for unit i is denoted x j,m i . Variable-blocks will be referred to through the corresponding matrix. T (resp. T 1 , ..., T p ) refers to a n × r T (resp. n × r 1 , ..., n × r p ) matrices of covariates. We assume that:

• The units, hence the rows of matrices Y, X 1 , ..., X p are independent multivariate normal vectors.

Model notations

The SEM we handle here is a restricted one, in that it contains only one structural equation, relating a dependent factor g, underlying a block Y of OV's, to p explanatory factors f 1 , ..., f p respectively underlying blocks X 1 , ..., X p of OV's (cf. fig. 1). The main assumptions of this model are the following:

• Factors f 1 , ..., f p are standard normal, i.e. ∀m ∈ 1, p , E(f m ) = 0 and V(f m ) = I n .
• In each block (e.g. X p ), the OV's (e.g. x m j , j ∈ 1, q p ) depend linearly on the block's factor (e.g. f m ) and a block of extra-covariates (e.g. T m ), conditional on which they are independent.

• Factor g is normal with zero-mean, and its expectation conditional on f 1 , ..., f p is a linear combination of them.

The SEM consists of p + 1 measurement equations and one structural equation. It is graphed on (cf. fig. 1).

Figure 1: Model with a dependent block and p explanatory blocks

Measurement equations

As formerly mentioned, each measurement equations relates the variables in a block X m (respectively Y ) to the block's factor f m (resp. g). This link may also involve covariates T m (resp. T ): each OV is expressed as a linear combination of the factor, the covariates and some noise. Hence the p + 1 measurement equations:

Y = T D + gb + ε Y ∀m ∈ 1, p , X m = T m D m + f m a m + ε m where D (resp. D m ) is a r T × q Y (resp. r m × q m ) parameter matrix, b (resp. a m ) a 1 × q Y (resp. 1 × q m )
parameter matrix, and ε Y (resp. ε m ) an n × q Y (resp. n × q m ) measurement-error matrix. We impose that the first column of T as well as of each T m matrix is equal to constant vector having all elements equal to 1. Thus, the first row of D and of each D m contains mean-parameters. As far as distributions are concerned, we assume that

• ε Y i ∼ N (0, ψ Y ) , where ψ Y = diag(σ 2 Y,j ) j∈ 1,q Y . • ∀m ∈ 1, p : ε m i ∼ N (0, ψ m )
, where ψ m = diag(σ 2 m,j ) j∈ 1,qm and that

• ε Y and ε m , ∀m ∈ 1, p are independent.
As to the factors, we assume that:

• ∀m ∈ 1, p : f m ∼ N (0, Id n ) with f 1 , ..., f m independent.

Structural equations

The structural equation we consider relates dependent factor g to explanatory factors f 1 , ..., f p (cf. fig. 1) through a linear model:

g = f 1 c 1 + • • • + f p c p + ε g
where ∀m ∈ 1, p , c m is a scalar parameter, and ε g ∈ R n is a disturbance vector.

We assume that

• ε g ∼ N (0, 1) • ε g is independent of ε Y and ε m , ∀m ∈ 1, p .
N.B. The unit-variance of disturbance ε g serves an identification purpose.

Hence we have the overall model:

   Y = T D + gb + ε Y ∀m ∈ 1, p , X m = T m D m + f m a m + ε m g = f 1 c 1 + • • • + f p c p + ε g (1)
where the set of parameters is

θ = {D, D 1 , ..., D p , b, a 1 , ..., a p , c 1 , c 2 , ψ Y , ψ 1 , ..., ψ p } such as θ is K-dimensional.
Thus, when all ψ matrices are diagonal, we have:

K = 2 + q Y (r T + 2) + p m=1 q m (r m + 2) (2)

A simplified model

But in order to avoid heavy formulas in the development of the algorithm, we shall use in the sequel, with no loss of generality, a simplified model involving p = 2 explanatory blocks X 1 and X 2 . The corresponding equation set, for a given unit i, reads:

       y i = t i D + g i b + ε i y x 1 i = t 1 i D 1 + f 1 i a 1 + ε i 1 x 2 i = t 2 i D 2 + f 2 i a 2 + ε i 2 g i = f 1 i c 1 + f 2 i c 2 + ε i g (3) Such as, θ = {D, D 1 , D 2 , b, a 1 , a 2 , c 1 , c 2 , σ 2 Y , σ 2 1 , σ 2 2 }.
Thus, in this case (cf. (2)), the dimension of θ is:

K = 5 + q Y (r T + 1) + 2 m=1 q m (r m + 1)

Estimation using the EM algorithm

In this work, likelihood maximization is carried out via an iterative EM algorithm [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF], section 4.7). Each iteration of the algorithm involves an Expectation (E)-step followed by a Maximization (M)-step. [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF] prove that the EM algorithm yields maximum likelihood estimates. Moreover, they proved that even if the starting point is one where the likelihood is not convex, if an instance of the algorithm converges, it will converge to a (local) maximum of the likelihood. Another major advantage of the EM algorithm is that it can be used to "estimate" missing data. Thus, if we consider LV's as missing data, the EM algorithm will prove a general technique to maximize the likelihood of statistical models with LV's, but also to estimate these LV's. In our SEM framework, LV's correspond to factors. Thus, we will be able to estimate the factors at unit-level. We shall present the algorithm on the simplified model (cf. section 2.4) with no loss of generality.

The EM algorithm

Let z = (y, x 1 , x 2 ) be the OV's and h = (g, f 1 , f 2 ) the LV's. The EM algorithm is based on the log-likelihood associated with the complete data (z, h).

The complete log-likelihood function

Let p(z, h; θ) denote the probability density of the complete data. The corresponding log-likelihood function is:

L(θ; z, h) = - 1 2 n i=1 {ln|ψ Y | + ln|ψ 1 | + ln|ψ 2 | + (y i -D t i -g i b) ψ -1 Y (y i -D t i -g i b) + (x 1 i -D 1 t 1 i -f 1 i a 1 ) ψ -1 1 (x 1 i -D 1 t 1 i -f 1 i a 1 ) + (x 2 i -D 2 t 2 i -f 2 i a 2 ) ψ -1 2 (x 2 i -D 2 t 2 i -f 2 i a 2 ) + (g i -c 1 f 1 i -c 2 f 2 i ) 2 + (f 1 i ) 2 + (f 2 i ) 2 } + λ
Where θ is the K-dimensional set of model parameters and λ a constant. However, because of the simplification made in the section 2.4, in our case θ

= {D, D 1 , D 2 , b, a 1 , a 2 , c 1 , c 2 , σ 2 Y , σ 2 1 , σ 2 2 }. Indeed, ψ Y = σ 2 Y Id q Y , ψ 1 = σ 2 1 Id q 1 and ψ 2 = σ 2 2 Id q 2 .

Estimation in the SEM

To maximize this function, in the framework of EM algorithm, we have to solve:

E h z [ ∂ ∂θ L(θ; z, h)] = 0 (4)
Foulley (2002) This demands that we know the derivatives of the log-likelihood function and the distribution p h i z i of h i conditional on z i for each observation i ∈ 1, n . Let us introduce the following notation:

p h i z i = N (M i =   m 1i m 2i m 3i   , Σ =   σ 11 σ 12 σ 13 σ 21 σ 22 σ 23 σ 31 σ 32 σ 33   ) g i = E h i z i [g i ] = m 1i ; γ i = E h i z i [g 2 i ] = (E h i z i [g i ]) 2 + V h i z i [g i ] = m 1 2 i + σ 11 f 1 i = E h i z i [f 1 i ] = m 2i ; φ 1 i = E h i z i [(f 1 i ) 2 ] = (E h i z i [f 1 i ]) 2 + V h i z i [f 1 i ] = m 2 2 i + σ 22 f 2 i = E h i z i [f 2 i ] = m 3i ; φ 2 i = E h i z i [(f 2 i ) 2 ] = (E h i z i [f 2 i ]) 2 + V h i z i [f 2 i ] = m 3 2 i + σ 33 For all ξ ∈ { g, f 1 , f 2 , γ, φ 1 , φ 2 }, we denote ξ = ( ξ i ) i=1,...,n ∈ R n .
The parameters of the gaussian distribution p h i z i are explicit and have the following form:

M i = Σ * 2 Σ * 3 -1 µ * and Σ = Σ * 1 -Σ * 2 Σ * 3 -1 Σ * 2 where: Σ * 1 =   (c 1 ) 2 + (c 2 ) 2 + 1 c 1 c 2 c 1 1 0 c 2 0 1   Σ * 2 =    ((c 1 ) 2 + (c 2 ) 2 + 1)b c 1 a 1 c 2 a 2 c 1 b a 1 0 (1,q 2 ) c 2 b 0 (1,q 1 ) a 2    Σ * 3 =    ((c 1 ) 2 + (c 2 ) 2 + 1)bb + Ψ Y c 1 ba 1 c 2 ba 2 c 1 a 1 b a 1 a 1 + Ψ 1 0 (q 1 ,q 2 ) c 2 a 2 b 0 (q 2 ,q 1 ) a 2 a 2 + Ψ 2    µ * i =   y i -D t i x 1 i -D 1 t 1 i x 2 i -D 2 t 2 i  
These results are demonstrated in AppendixB. Expressions of the first-order derivatives of L with respect to θ are also established in AppendixB and written in the following forms with m ∈ {1, 2}:

                                                 ∂ ∂D L(z, h) = n i=1 ψ -1 Y (y i -D t i -g i b)t i ∂ ∂D m L(z, h) = n i=1 ψ -1 m (x m i -D m t m i -f m i a m )t m i ∂ ∂b L(z, h) = n i=1 g i ψ -1 Y (y i -D t i -g i b) ∂ ∂a m L(z, h) = n i=1 f m i ψ -1 m (x m i -D m t m i -f m i a m ) ∂ ∂c m L(z, h) = n i=1 f m i (g i -c 2 f 2 i -c 1 f 1 i ) ∂ ∂σ 2 Y L(z, h) = n q Y σ -2 Y -σ -4 Y n i=1 ||y i -D t i -g i b|| 2 ∂ ∂σ 2 m L(z, h) = n q m σ -2 m -σ -4 m n i=1 ||x m i -D m t m i -f m i a m || 2 (5)
So, here formula (4) develops into:

                                                         n i=1 (y i -D t i -g i b)t i = n i=1 (x m i -D m t m i -f m i a m )t m i = n i=1 g i y i -g i D t i -γ i b = n i=1 f m i x m i -f m i D m t m i -φ m i a m = n i=1 σ 12 + f 1 i g i -c 2 σ 23 -c 2 f 1 i f 2 i -φ 1 i c 1 = n i=1 σ 31 + f 2 i g i -c 2 φ 2 i -c 1 σ 32 -c 1 f 1 i f 2 i = nq Y σ -2 Y -σ -4 Y n i=1 ||y i -D t i || 2 + ||b|| 2 γ i -2(y i -D t i ) g i b = nq m σ -2 m -σ -4 m n i=1 ||x m i -D m t m i || 2 + ||a m || 2 φ m i -2(x m i -D m t m i ) f m i a m = (6)
System of equations ( 6) is easy to solve and the obtained solutions will be given in the next section.

Results

The explicit solution of the system (6) (and also of ( 4)) is the following:

b = ( gy -yt )(tt ) -1 gt γ -gt (tt ) -1 gt a m = f m x m -x m t m (t m t m ) -1 f m t m φ m -f m t m (t m t m ) -1 f m t m c 1 = (σ 12 + f 1 g) φ 2 -(σ 13 + f 2 g)(σ 23 + f 1 f 2 ) φ 1 φ 2 -(σ 23 + f 1 f 2 ) 2 c 2 = (σ 13 + f 2 g) φ 1 -(σ 12 + f 1 g)(σ 23 + f 1 f 2 ) φ 1 φ 2 -(σ 23 + f 1 f 2 ) 2 D = (yt -b gt )(tt ) -1 D m = (x m t m -a m f m t m )(t m t m ) -1 σ 2 Y = 1 nq Y n i=1 {||y i -D t i || 2 + || b|| 2 γ i -2(y i -D t i ) b g i } σ 2 m = 1 nq m n i=1 {||x m i -D m t m i || 2 + || a m || 2 φ m i -2(x m i -D m t m i ) a m f m i } (7)

The algorithm

To estimate parameters in θ, we propose the following EM-algorithm. We denote [t] the t ieth -iteration of the algorithm.

1. Initialization 1 = choice of the initial parameter values θ [0] .

1 In the initialization step, ∀m ∈ {1, p} we propose to obtain D m[0] by multiple linear regression between X m and T m . Then, to initialize the others, we compute each approximated factor f m [0] and g [0] as first principal component of a PCA of X m -T m D m[0] and Y -T D [0] . Thus, we initialize a m , σ 2 m (resp. b, σ 2 y ) by multiple linear regression between X m -T m D m[0] and f m [0] (resp. between Y -T D [0] and g [0] ). Finally, each c m[0] can be obtained by multiple linear regression between g [0] and

p m=1 f m [0]
. In practice we use the functions lm() and PCA() derived from the package FactoMineR [START_REF] Husson | FactoMineR: An R Package for Multivariate Analysis[END_REF].

Current iteration

t ≥ 1, until stopping condition is met: (a) E-step: with θ [t-1] , i. Calculate explicitly distribution p h i z i for each i ∈ 1, n . ii. Estimate the factor-values g [t] , f m [t]
, m ∈ {1, 2}.

iii. Calculate γ [t] and

φ m [t] , m ∈ {1, 2}. (b) M-step:
i. Update θ to θ [t] by injecting g [t] , γ [t] and f m [t] ,

φ m [t]
, m ∈ {1, 2} into the formulas in (7). 3. We used the following stopping condition with the smallest possible:

K k=1 |θ * [t+1] [k] -θ * [t] [k]| θ * [t+1] [k] < (8)
where θ * is the K-dimensional vector containing the scalar values in all parameters in θ.

Numerical results on simulated data

Data generation

We consider n = 400 units and q Y = q 1 = q 2 = 40. Therefore, the 120 OV's Y, X 1 , X 2 are simulated so as to be structured respectively around three factors g, f 1 , f 2 . Factors f 1 and f 2 are explanatory of g. Besides, we consider r T = r 1 = r 2 = 2 i.e 2 covariates are simulated for each covariate matrix T , T 1 and T 2 . The data is simulated as follows.

Choice of θ:

(a) D = D 1 = D 2 a) matrices filled in row-wise with the ordered integer sequence ranging from 1 to 80 (indeed:

r T * q Y = r 1 * q 1 = r 2 * q 2 = 2 * 40). (b) b = a 1 = a 2 = ordered integer sequence ranging from 1 to 40. (c) c 1 = c 2 = 1 (d) σ 2 Y = σ 2 1 = σ 2 2 = 1 2. Simulation of factors g, f 1 , f 2
(a) Simulate vectors f 1 and f 2 of n = 400 normally distributed random numbers with mean 0 and variance 1 (abbreviated ∀m, ∈ {1, 2} f m ∼ N (0, Id 400 )).

(b) We simulate ε g according to distribution ε g ∼ N (0, Id 400 ).

(c) We then calculate g as

g = f 1 c 1 + f 2 c 2 + ε g 3. Simulation of noises ε Y , ε 1 , ε 2 Each element of matrix ε Y , (respectively ε 1 , ε 2 ) is simulated independently from distribution N (0, σ 2 Y = 1) (re- spectively, σ 2 1 = 1, σ 2 2 = 1). 4. Simulation of covariate matrices T , T 1 , T 2
Each element of matrices T , T 1 , T 2 is simulated according to the standard normal distribution. 5. Construction of Y , X 1 , X 2 Y , X 1 , X 2 are eventually calculated through formulas in the model (1).

This simulation scheme was performed 100 times, each time yielding a set of simulated data matrices (Y, X 1 , X 2 ). Then for each simulated data, we ran an estimation routine with a threshold value ε = 10 -2 , yielding the average results presented in section 4.2. Thus from 400 * 120 = 48000 scalar elements of data, we will estimate 3 * n = 1200 scalar elements of factors plus K = 5 + 3 * 40(2 + 1) = 365 scalar parameters, i.e: 1565 scalars.

Results

Convergence was observed in almost all cases in less than five iterations. We assess the quality of the estimations as follows.

• On the one hand, we calculate the absolute relative deviation between each simulated scalar parameter in θ * and its estimation, and then average these deviations over the 100 simulations. We then produce a box-plot of the average absolute relative deviations (cf. fig. 2). This makes the interpretation easier, since we only need to look at the boxplot's values and check that they are positive (because of the absolute value) and close to 0.

• On the other hand, to assess the quality of the factor estimations, we compute the 300 values of square correlations between the simulated concatenated factors (g, f 1 , f 2 ) (respectively) and the corresponding estimations (( g, f 1 , f 2 )). Once again, we produce a box-plot of these correlations (cf. fig. 3) and check that it indicates values close to 1.

Figures 2 and3 show clearly that the estimations are very close to the actual quantities. Indeed, on figure 2, the median of average absolute relative deviations is 0.018, the first quartile is 0.015 and the third quartile is 0.023. On figure 3, the median of square correlations is 0.998, the first quartile is 0.997 and the third quartile is 0.999. So, factor g (respectively f 1 and f 2 ) turn out to be drawn towards the principal direction underlying the bundles made up by observed variables Y (respectively X 1 and X 2 ). Now, we may legitimately wonder how the quality of estimations could be affected by the number of observations and the number of OV's in each block. In the following section we give a sensitivity analysis performed to investigate this.

Sensitivity analysis of estimations

We performed a sensibility analysis on the simulated data presented in section 4.1. The purpose was to study the influence of the block-sizes (n, q Y , q 1 , q 2 ) on the quality of estimation, both of the parameters and the factors. To simplify the analysis, we imposed q Y = q 1 = q 2 = q and varied n and q separately, i.e. studied the cases n = 50, 100, 200, 400 with q = 40 and q = 5, 10, 20, 40 with n = 400. Each case was simulated 100 times. Therefore, we simulated 800 data-sets.

Sensitivity with respect to the number n of observations

In this section, we study the evolution with n of the average estimation of structural coefficients c 1 and c 2 and parameter σ 2

Y with respect to their actual values, all equal to 1, and that of the correlations between factors and their estimates. The number of OV's is fixed to q = 40 in each block. Figures 4,5 and 6 graph these evolutions (average value of estimate in plain line), including a 95% confidence-interval about each average estimate (dotted line). These figures show that the biases and the standard deviations are, as expected, more important for little values of n, but also that the quality of estimation is already quite good for n = 50. As for the factors, figure 7 shows that their correlations increase and get close to one as n increases, with a dispersion decreasing to 0. However, even for n = 50, the correlations are mostly above 0.95, indicating that the factors are correctly reconstructed.

Sensitivity with respect to the number q of OV's in each block

Likewise, we study the evolution of the average estimates of c 1 , c 2 , σ 2 Y and the correlation between factors and their estimates for different values of q, with n fixed to 400. We observe that, unsurprisingly, the biases and the standard deviations decrease as q increases (cf. figures 8, 9 and 10). We observe that they stabilize even faster with q than with n, particularly σ 2 Y . Indeed, from q = 10 on, the confidence interval is narrow enough. As for the factors, figure 11 shows that their correlations are already very close to 1 for q = 5, with a very small variance, and keep increasing with q. To sum things up, the sample size n proved to have more impact on the quality of parameter estimation and factor reconstruction than the number of OV's. Now, the quality of factor reconstruction remains high for rather small values of n or q. We advise to use a minimal sample size of n = 100 to obtain really stable structural coefficients. Above this threshold, n has but little impact on the biases and standard deviations of estimations.

An application to environmental data

Data presentation

We apply our model to the data-set genus provided in the R-package SCGLR by [START_REF] Mortier | SCGLR-An R Package for Supervised Component Generalized Linear Regression[END_REF]. Data-set genus was built from the CoFor-Change database. It gives the abundances of 27 common tree genera present in the tropical moist forest of the Congo-Basin, and the measurements of 40 geo-referenced environmental variables, for n = 1000 inventory plots (observations). Some of the geo-referenced environmental variables describe 16 physical factors pertaining to topography, geology and rainfall description. The remaining variables characterize vegetation through the enhanced vegetation index (EVI) measured on 16 dates. In this section, we aim at modeling the tree abundances from the other variables, while reducing the dimension of data. The dependent block of variables Y therefore consists of the q Y = 27 tree species counts divided by the plot-surface. A PCA of the geo-referenced environmental variables and the photosynthetic activity variables confirms that EVI measures are clearly separated from the other variables (cf. Fig. 12). Indeed, Fig. 12 shows two variable-bundles with almost orthogonal central directions. This justifies using our model (cf. section 5.2) with p = 2 explanatory groups, one of them (X 1 ) gathering q 1 = 16 rainfall measures and location variables (longitude, latitude and altitude), and the second one (X 2 ), the q 2 = 23 EVI measures. Besides, in view of the importance of the geological substrate on the spatial distribution of tree species in the Congo Basin, showed by [START_REF] Fayolle | Geological Substrates Shape Tree Species and Trait Distributions in African Moist Forests[END_REF], we chose to put nominal variable geology in a block T . This block therefore contains constant 1 plus all the indicator variables of geology but one, which will be the reference value. Geology having 5 levels, T has thus 5 columns.

Model with geologic covariates 5.2.1. Model specification

Here is the model used with the variable-blocks designed in section 5.1.:

       Y = T D + gb + ε Y X 1 = 1 n d 1 + f 1 a 1 + ε 1 X 2 = 1 n d 2 + f 2 a 2 + ε 2 g = f 1 c 1 + f 2 c 2 + ε g
Where n = 1000, q Y = 27, q 1 = 16, q 2 = 23 and r T = 5. The first row of D is a parameter vector that contains the means of the Y 's noted D[1, ] in Table 1, and the other rows, the overall effects of the geological substrates with respect to the reference one. Indeed, the next section presents the model's parameter-estimations where, in Table 1, each row r of D is noted D[r, ].

Results

With a threshold value ε = 10 -3 , convergence was reached after 58 iterations. Some parameter-estimations are presented in Tables 1, 2 and3 It can be seen in Tables 1 and6 that for certain species, the geologic substrate seems to be of great importance (e.g. for gen1, gen5, gen7, gen9, gen12, gen16, gen21, gen25, gen26, gen27), whereas for others, it only has a small impact on the abundances (e.g. for gen2, gen6, gen8, gen10, gen18, gen20, gen23). Moreover, Table 1 shows that the correlations between g and Y are high in absolute value only for few variables : gen2, gen23, gen24 and gen25. Therefore, only these are well accounted for by our model. Although we have carried out the analysis with variables gen2, gen3, gen8, gen10, gen11, gen15, gen17, gen23, gen24 and gen25, the results are practically the same when we take all variables. In table 2, the correlations between f 1 and variables pluvio_1 to pluvio_12 of X 1 show two rainfall regimes. Indeed, pluvio_1 corresponds to january, pluvio_2, to february, ... pluvio_12 to december. The Central African Republic has a tropical climate : the dry season ranges from November to April and the rainy season from June to September. Figure 13 shows that f 1 is positively correlated to the rainfalls of the rainy season and negatively to those of the dry one.

Model without covariate

What if we omit the geologic substrate as covariate?

Model specification

We now consider the model without covariate geology, i.e. every T is reduced to

1 n :        Y = 1 n d + gb + ε Y X 1 = 1 n d 1 + f 1 a 1 + ε 1 X 2 = 1 n d 2 + f 2 a 2 + ε 2 g = f 1 c 1 + f 2 c 2 + ε g
Where, there always are n = 1000, q Y = 27, q 1 = 16, q 2 = 23. Also, d (resp. d 1 and d 2 ) are parameters vectors that contains the means of the Y 's (resp. X 1 and X 2 ).

Results

With a threshold value ε = 10 -3 , convergence was reached after 49 iterations. Some parameter-estimations are presented in Table 4 and5 Table 3 and 5 show that, the geological effect is being considered or removed, the geographic factor (position and rainfalls) keeps a much greater effect than the EVI's. The estimations of σ 2 1 , σ 2 2 , σ 2 Y don't change significantly. For both, the rainfall regimes are identically identified. However, table 2 shows the great impact of geologic substrate on the abundance of species gen1, gen5, gen7, gen9, gen12, gen16, gen21, gen25, gen26 and gen27. Therefore, the presence of covariate geology in the model is relevant.

Assessing the model quality through re-sampling

To assess the stability of results and thus, validate the models (with covariate), we use a 5-fold re-sampling technique: 5 separate samples are randomly extracted from the complete genus data, thus, their size is kept equal to 200 units. For each of them, we obtain estimated parameters and factors. Then, for each sample, we compute an average Mean Square Error (MSE) and an average correlation between the parameter-estimates obtained on the sample and those obtained on the complete data. Finally, on each sample, we calculate an average MSE and correlation between the factorestimates obtained on the sample and the corresponding ones obtained on the complete data for the units belonging to the sample.

Model with geologic covariate

Figure 14 (resp. 15) shows the average MSE (resp. the correlation) between estimated parameters on 5 data samples θ * s∈ 1,5 and estimated parameters on the complete data θ * . More precisely, for these average MSE (respectively correlation), the median is 3.85 * 10 -3 (resp. 0.99), the first quartile is 1.95 * 10 -3 (resp. 0.99) and the third quartile is 6.17 * 10 -3 (resp. 0.99). These values are close to 0 (resp. 1). So, we can be rather confident in the estimates of parameters obtained in the last section.

Conclusion

The maximum-likelihood estimation method is known to be a stringent method of estimation having nice properties. In the context of estimation methods of SEM, the LISREL approach is based on likelihood maximization, contrary to PLS, THEME, and other component-based methods. However, LISREL only focusses on the variance-covariance structure of the data and does not allow to estimate the LV's, contrary to PLS and THEME. To estimate them, we proposed to carry out likelihood maximization through the EM algorithm. This approach assumes that VLs are factors, which is constraining than assuming they are components. Therefore this approach combines the stringency of likelihood maximization with the possibility to estimate the LV's. This presented this new approach and performed sensitivity analysis to show its performances. Then, an application on environmental data was made, which shows how to use this method. Along with re-sampling for validation purposes.

Acknowledgments

The forest inventories were funded by the French Agency for Development (AFD) through the PARPAF project. We wish to thank S. Chong (TCA), A. Banos (SCAF), and the '"Ministère des Eaux, Forêts, Chasse et Pêche" of the Central African Republic for authorizing access to the inventory data, and the field teams who drew up these inventories. This study is part of the ErA Net BiodivERsA CoForChange project, funded by the National Research Agency (ANR) and the Natural Environment Research Council (NERC), involving 16 European, African and international partners and a number of timber companies (see the list on the website, http://www.coforchange.eu).

Appendix A. Calculation of the complete data log-likelihood function L Proof. In our case p = 2, ψ Y = σ 2 Y Id q Y , ψ 1 = σ 2 1 Id q 1 and ψ 2 = σ 2 2 Id q 2 , and for observation i, the model is formulated as follows:

       y i = t i D + g i b + ε i y x 1 i = t 1 i D 1 + f 1 i a 1 + ε i 1 x 2 i = t 2 i D 2 + f 2 i a 2 + ε i 2 g i = f 1 i c 1 + f 2 i c 2 + ε i g We have, p(z i , h i ; θ) = p(y i , x 1 i , x 2 i , g i , f 1 i , f 2 i ; θ) = p(y i , x 1 i , x 2 i |g i , f 1 i , f 2 i ; θ)p(g i , f 1 i , f 2 i ; θ) = p(y i , x 1 i , x 2 i |g i , f 1 i , f 2 i ; θ)p(g i |f 1 i , f 2 i ; θ)p(f 1 i , f 2 i ; θ) = p(y i , x 1 i , x 2 i |g i , f 1 i , f 2 i ; θ)p(g i |f 1 i , f 2 i ; θ)p(f 1 i ; θ)p(f 2 i ; θ) = p(y i , x 1 i , x 2 i |g i , f 1 i , f 2 i ; θ)p(g i |f 1 i , f 2 i ; θ)p(f 1 i )p(f 2 i ) = p(x 1 i , x 2 i |y i , g i , f 1 i , f 2 i ; θ)p(y i |g i , f 1 i , f 2 i ; θ)p(g i |f 1 i , f 2 i ; θ)p(f 1 i )p(f 2 i ) = p(x 1 i , x 2 i |f 1 i , f 2 i ; θ)p(y i |g i ; θ)p(g i |f 1 i , f 2 i ; θ)p(f 1 i )p(f 2 i ) = p(x 1 i |x 2 i , f 1 i , f 2 i ; θ)p(x 2 i |f 1 i , f 2 i ; θ)p(y i |g i ; θ)p(g i |f 1 i , f 2 i ; θ)p(f 1 i )p(f 2 i ) = p(x 1 i |f 1 i ; θ)p(x 2 i |f 2 i ; θ)p(y i |g i ; θ)p(g i |f 1 i , f 2 i ; θ)p(f 1 i )p(f 2 i ) Where θ = {D, D 1 , D 2 , b, a 1 , a 2 , c 1 , c 2 , ψ Y , ψ 1 , ψ 2 } is the set of model param- eters. Therefore, L(θ; z i , h i ) = L(θ; x 1 i |f 1 i ) + L(θ; x 2 i |f 2 i ) + L(θ; y i |g i ) + L(θ; g i |f 1 i , f 2 i ) + L(f 1 i ) + L(f 2 i )
Because of the model and the normal distribution properties we obtain:

x m i |f m i ∼ N (t m i D m + f m i a m , ψ X m ) y i |g i ∼ N (t i D + g i b , ψ Y ) g i |f 1 i , f 2 i ∼ N (f 1 i c 1 + f 2 i c 2 , 1) f m i ∼ N (0, 1)
Then, we obtain the complete data log-likelihood function:

L(θ; z, h) = - 1 2 n i=1 {ln|ψ Y | + ln|ψ 1 | + ln|ψ 2 | + (y i -D t i -g i b) ψ -1 Y (y i -D t i -g i b) + (x 1 i -D 1 t 1 i -f 1 i a 1 ) ψ -1 1 (x 1 i -D 1 t 1 i -f 1 i a 1 ) + (x 2 i -D 2 t 2 i -f 2 i a 2 ) ψ -1 2 (x 2 i -D 2 t 2 i -f 2 i a 2 ) + (g i -c 1 f 1 i -c 2 f 2 i ) 2 + (f 1 i ) 2 + (f 2 i ) 2 } + λ
Where λ a constant. Also, the set of model parameters

θ = {D, D 1 , D 2 , b, a 1 , a 2 , c 1 , c 2 , ψ Y , ψ 1 , ψ 2 } in our case corresponds to θ = {D, D 1 , D 2 , b, a 1 , a 2 , c 1 , c 2 , σ 2 Y , σ 2 1 , σ 2
2 } because of the simplification made in the section 2.3, in our case. Indeed, ψ Y = σ 2 Y Id q Y , ψ 1 = σ 2 1 Id q 1 and ψ 2 = σ 2 2 Id q 2 . Therefore, we can also write the complete data log-likelihood function:

L(θ; z, h) = - 1 2 n i=1 {q Y ln(σ 2 Y ) + q 1 ln(σ 2 1 ) + q 2 ln(σ 2 2 ) + σ -2 Y (y i -D t i -g i b) (y i -D t i -g i b) + σ -2 1 (x 1 i -D 1 t 1 i -f 1 i a 1 ) (x 1 i -D 1 t 1 i -f 1 i a 1 ) + σ -2 2 (x 2 i -D 2 t 2 i -f 2 i a 2 ) (x 2 i -D 2 t 2 i -f 2 i a 2 ) + (g i -c 1 f 1 i -c 2 f 2 i ) 2 + (f 1 i ) 2 + (f 2 i ) 2 } + λ Appendix B. Demonstration of the normality of the distribution of h i |z i Proof. In our case p = 2, ψ Y = σ 2 Y Id q Y , ψ 1 = σ 2 1
Id q 1 and ψ 2 = σ 2 2 Id q 2 , and for observation i, the model is formulated as follows:

       y i = t i D + g i b + ε i y x 1 i = t 1 i D 1 + f 1 i a 1 + ε i 1 x 2 i = t 2 i D 2 + f 2 i a 2 + ε i 2 g i = f 1 i c 1 + f 2 i c 2 + ε i g
To prove the normality of the distribution of h i |z i presented in section 3.1.2., we use the classical result 2 about the conditioning of normally distributed variables. Before using this result, we calculate the joint distribution of (g

i , f 1 i , f 2 i , y i , x 1 i , x 2 i ).
We know that, for observation i,

y i ∼ N (D t i , b((c 1 ) 2 + (c 2 ) 2 + 1)b + Ψ Y ) x m i ∼ N (D m t m i , a m a m + Ψ m ) g i ∼ N (0, (c 1 ) 2 + (c 2 ) 2 + 1) f m i ∼ N (0, 1)
Then, after compute the required covariances we obtain,

(g i , f 1 i , f 2 i ) ∼ N (   0 0 0   ,   (c 1 ) 2 + (c 2 ) 2 + 1 c 1 c 2 c 1 1 0 c 2 0 1   )
And,

(y i , x 1 i , x 2 i ) ∼ N (   D t i D 1 t 1 i D 2 t 2 i   ,    ((c 1 ) 2 + (c 2 ) 2 + 1)bb + Ψ Y c 1 ba 1 c 2 ba 2 c 1 a 1 b a 1 a 1 + Ψ 1 0 (q 1 ,q 2 ) c 2 a 2 b 0 (q 2 ,q 1 ) a 2 a 2 + Ψ 2    ) 2 If two variables X 1 and X 2 are normally distributed such that, X 1 X 2 ∼ N (µ = µ 1 µ 2 , Σ = Σ 11 Σ 12 Σ 21 Σ 22 )
where, µ 1 (r × 1), µ 2 (s × 1), Σ 11 (r × r), Σ 12 (r × s), Σ 21 (s × r) and Σ 22 (s × s); then,

(X 1 |X 2 = x 2 ) ∼ N (M = µ 1 + Σ 12 Σ 22 -1 (x 2 -µ 2 ), φ = Σ 11 -Σ 12 Σ 22 -1 Σ 21 ) (9)
Then, after compute the required covariances we obtain the joint distribution, (g

i , f 1 i , f 2 i , y i , x 1 i , x 2 i ) ∼ N (M * i , Σ * ) such as, M * i =     0 (3,1) D t i D 1 t 1 i D 2 t 2 i     and Σ * = Σ * 1 Σ * 2 Σ * 2 Σ *

3

. Where,

Σ * 1 =   (c 1 ) 2 + (c 2 ) 2 + 1 c 1 c 2 c 1 1 0 c 2 0 1   Σ * 2 =    ((c 1 ) 2 + (c 2 ) 2 + 1)b c 1 a 1 c 2 a 2 c 1 b a 1 0 (1,q 2 ) c 2 b 0 (1,q 1 ) a 2    Σ * 3 =    ((c 1 ) 2 + (c 2 ) 2 + 1)bb + Ψ Y c 1 ba 1 c 2 ba 2 c 1 a 1 b a 1 a 1 + Ψ 1 0 (q 1 ,q 2 ) c 2 a 2 b 0 (q 2 ,q 1 ) a 2 a 2 + Ψ 2   
Finally, we use result ( 9) and obtain the distribution,

h i |z i ∼ N (M i , Σ)
where,

M i = Σ * 2 Σ * 3 -1 µ * i and Σ = Σ * 1 -Σ * 2 Σ * 3 -1 Σ * 2 , such that µ * i =   y i -D t i x 1 i -D 1 t 1 i x 2 i -D 2 t 2 i  
Appendix C. Calculation of the first-order derivatives of L Proof. We search the first-order derivatives of the complete data log-likelihood function:

L(θ; z, h) = - 1 2 n i=1 {ln|ψ Y | + ln|ψ 1 | + ln|ψ 2 | + (y i -D t i -g i b) ψ -1 Y (y i -D t i -g i b) + (x 1 i -D 1 t 1 i -f 1 i a 1 ) ψ -1 1 (x 1 i -D 1 t 1 i -f 1 i a 1 ) + (x 2 i -D 2 t 2 i -f 2 i a 2 ) ψ -1 2 (x 2 i -D 2 t 2 i -f 2 i a 2 ) + (g i -c 1 f 1 i -c 2 f 2 i ) 2 + (f 1 i ) 2 + (f 2 i ) 2 } + λ
Where λ constant, θ = {D, D 1 , D 2 , b, a 1 , a 2 , c 1 , c 2 , ψ Y , ψ 1 , ψ 2 }, ψ Y = σ 2 Y Id q Y , ψ 1 = σ 2 1 Id q 1 and ψ 2 = σ 2 2 Id q 2 . Therefore, there are matrix-parameters (D, D 1 , D 2 ), vector-parameters (b, a 1 , a 2 ) and scalar parameters (c 1 , c 2 , σ 2 Y , σ 2 1 , σ 2 2 ). Then, L is a sum of three types of functions: the logarithm, the square function and a quadratic form function (w -Xβ) Γ(w -Xβ), where Γ is symmetric and w (q × 1), X (q × m), β (m × 1) and Γ (q × q). The first-order derivatives of the logarithm function and the square function are in our case trivial. The first-order derivative of (w -Xβ) Γ(w -Xβ) by X is less trivial but necessary. Let us start by making explicit the first-order derivative of (w -Xβ) Γ(w -Xβ) with respect to X. Likewise, we establish that :

∂ ∂D L(z, h) = n i=1 ψ -1 Y (y i -D t i -g i b)t i
Similar reasoning can be applied to D m and allows to obtain the second row of (5). Concerning the third and the fourth row of (5), we use the classical result :

∂ ∂β [(w -Xβ) Γ(w -Xβ)] = -2X Γ(w -Xβ)
Eventually, the fifth, the sixth and the eighth rows of (5) are obtained in a trivial way.

Parameter

Figure 2 :

 2 Figure 2: Box plot of the average absolute relative deviations between simulated parameters and their estimations.

Figure 3 :

 3 Figure 3: Box plot of the correlations between the simulated factors and their estimations

Figure 4 :

 4 Figure 4: Average estimate of c 1 and 95% confidence interval as a function of n.

Figure 5 :

 5 Figure 5: Average estimate of c 2 and 95% confidence interval as a function of n.

Figure 6 :

 6 Figure 6: Average estimate of σ 2 Y and 95% confidence interval as a function of n.

Figure 7 :

 7 Figure 7: Box plots of the correlations between simulated factors and their estimations for various values of n.

Figure 8 :

 8 Figure 8: Evolution of the average of estimated parameters c 1 according to various n values.

Figure 9 :

 9 Figure 9: Evolution of the average of estimated parameters c 2 according to various n values.

Figure 10 :

 10 Figure 10: Evolution of the average of estimated parameters σ 2 Y (at right) according to various n values.

Figure 11 :

 11 Figure 11: Box plots of the correlations between simulated factors and their estimations according to various n values.

Figure 12 :

 12 Figure 12: Correlation-scatterplot yielded by the PCA of the X 1 and X 2 geo-referenced environmental variables.

Figure 13 :

 13 Figure 13: Correlations between f 1 and the monthly variables of X 1 : two rainfall regimes.

Figure 14 :

 14 Figure 14: Box plot of the average MSE's between the parameter-estimates obtained on the 5 genus data sub-samples and those obtained on the complete data.

Figure 15 :

 15 Figure 15: Box plot of the average correlations between the parameter-estimates obtained on the 5 genus data sub-samples and those obtained on the complete data.

Figure 16 and

 16 Figure16and Figure17respectively give the box-plot of the factors' average MSE and correlation for each of the 5 samples. More precisely, for these average MSE's (respectively correlations), the median is 1.15 * 10 -2 (resp. 0.98), the first quartile is 7.44 * 10 -3 (resp. 0.98) and the third quartile is 3.53 * 10 -2 (resp. 0.99). These values are close enough to 0 (resp. 1) to allow us to be confident in the estimates obtained on the complete data.

Figure 16 :

 16 Figure 16: Box-plot of the average MSE of factor-estimates.

Figure 17 :

 17 Figure 17: Box-plot of the average correlation of factor-estimates.

  d X [(w -Xβ) Γ(w -Xβ)] = (w -Xβ) Γ(-dXβ) + (-dXβ) Γ(w -Xβ) = -2(w -Xβ) Γ(dXβ) = tr[-2(w -Xβ) Γ(dXβ)] = tr[-2β(w -Xβ) ΓdX] =< -2β(w -Xβ) Γ|dX > Therefore, d dX [(w -Xβ) Γ(w -Xβ)] = (-2β(w -Xβ) Γ) = -2(β(w -Xβ) Γ) = -2Γ(w -Xβ)β

Table 1 :

 1 Application to the genus data with geologic covariate : estimations of parameters D and b , and correlations between g and the variables Y

	.

Table 2 :

 2 Application to the genus data with geologic covariate : estimations of parameters d 1 and a 1 , and correlations between f 1 and the variables X 1

Table 3 :

 3

Application to genus data with geologic covariate: scalar parameter-estimations

Table 4 :

 4 Application to the genus data without covariate : estimations of parameters D and b , and correlations between g and the variables Y

	. For

Table 5 :

 5 Application to the genus data without covariate : scalar estimations of parameters

Table 7 :

 7 Application to the genus data with geologic covariate : estimations of parameters d 2 and a 2 , and correlations between f 2 and the variables X 2

			-estimations
	Variables	d 2	a 2	Correlations with f 2
	evi_1	15.51	0.63	0.65
	evi_2	13.47	0.59	0.6
	evi_3	14.83	0.51	0.52
	evi_4	14.67	0.58	0.6
	evi_5	16.44	0.56	0.57
	evi_6	18.74	0.51	0.52
	evi_7	18.44	0.75	0.76
	evi_8	20.59	0.8	0.82
	evi_9	21.83	0.76	0.78
	evi_10	19.19	0.74	0.76
	evi_11	18.22	0.67	0.69
	evi_12	15.92	0.61	0.63
	evi_13	15.4	0.58	0.6
	evi_14	13.51	0.7	0.72
	evi_15	14.57	0.69	0.71
	evi_16	14.95	0.76	0.78
	evi_17	16.09	0.73	0.75
	evi_18	15.95	0.77	0.79
	evi_19	17.12	0.73	0.75
	evi_20	15.02	0.75	0.77
	evi_21	15.87	0.75	0.77
	evi_22	14.21	0.71	0.73
	evi_23	15.26	0.68	0.69

Table 8 :

 8 Application to the genus data without covariate : estimations of parameters d 1 and a 1 , and correlations between f 1 and the variables X 1

		Parameter-estimations
	Variables	d 2	a 2	Correlations with f 2
	evi_1	15.51	0.63	0.65
	evi_2	13.47	0.59	0.6
	evi_3	14.83	0.51	0.52
	evi_4	14.67	0.58	0.6
	evi_5	16.44	0.56	0.57
	evi_6	18.74	0.5	0.52
	evi_7	18.44	0.74	0.76
	evi_8	20.59	0.8	0.82
	evi_9	21.83	0.76	0.78
	evi_10	19.19	0.74	0.76
	evi_11	18.22	0.67	0.69
	evi_12	15.92	0.61	0.63
	evi_13	15.4	0.58	0.6
	evi_14	13.51	0.7	0.72
	evi_15	14.57	0.69	0.71
	evi_16	14.95	0.76	0.78
	evi_17	16.09	0.73	0.75
	evi_18	15.95	0.77	0.79
	evi_19	17.12	0.73	0.75
	evi_20	15.02	0.75	0.77
	evi_21	15.87	0.75	0.77
	evi_22	14.21	0.71	0.73
	evi_23	15.26	0.67	0.69

Table 9 :

 9 Application to the genus data without geologic covariate : estimations of parameters d 2 and a 2 , and correlations between f 2 and the variables X 2

Appendix D. Tables of section 5.

Variables


Prométhée Spathis
email: promethee.spathis@lip6.fr

Salah-Eddine Belouanas
email: salah-eddine.belouanas@lip6.fr

Kim-Loan Thai
email: kim.thai@lip6.fr

Leveraging Replication in Content-Centric Networks

This paper focuses on the communication functionalities of content-centric networks. More especially, we take advantage of the widespread availability of content as today data objects may be replicated or moved at several locations to minimize content access time or to provide resiliency. We propose a scalable interdomain routing architecture designed to meet the needs of data-intensive applications. Our protocol complements CCN's routing framework with the ability to discover routes that allow users to retrieve replicated content from the most efficient location. The distinctive feature of our approach lies in the use of the caching capabilities of content routers to store and maintain the routes discovered towards the closest replicas. Our simulation results show an increase in content access efficiency by minimizing retrieval time and by maximizing the bandwidth demand.

I. INTRODUCTION

In the last years, there has been a revival of interest in content-centric networking that have lead to the proposal of several architectures following popular designs as [START_REF] Cheriton | An architecture for content routing support in the internet[END_REF], [START_REF] Caesar | Rofl: routing on flat labels[END_REF], [START_REF] Koponen | A data-oriented (and beyond) network architecture[END_REF]. Motivated by the ability to access content in a locationindependent manner, the common underlying paradigm to these architectures calls for a shift in the way communication takes place among network entities. The emphasis is given to the content rather than the location. As such, the principles behind the design of content-centric architectures fit within the broader framework of Information-Centric Networking (ICN) [START_REF] Xylomenos | A survey of information-centric networking research[END_REF] which is also motivated by the challenges arising from the current end-to-end host-centric design of the Internet.

The most relevant ICN architectures include the seminal papers of Gritterand and Cheriton [START_REF] Cheriton | An architecture for content routing support in the internet[END_REF] and Carzaniga et al. [START_REF] Carzaniga | Forwarding in a content-based network[END_REF], NSF projects like NDN [START_REF] Zhang | Named data networking (ndn) project[END_REF] and MobilityFirst [START_REF]Nsf mobility first project[END_REF], the European projects PSIRP [START_REF]Fp7 psirp project[END_REF], PURSUIT [START_REF]Fp7 pursuit project[END_REF], 4WARD [START_REF]Fp7 4ward project[END_REF], COMET [START_REF]Fp7 comet project[END_REF], and SAIL [START_REF]Fp7 sail project[END_REF], as well as the French funded project ANR Connect [START_REF]Anr connect project[END_REF] which adopts the NDN architecture. ICN architectures address research issues such as naming, innetwork caching, security and mobility [START_REF] Carofiglio | From content delivery today to information centric networking[END_REF], [START_REF] Diallo | A contentbased publish/subscribe framework for large-scale content delivery[END_REF], [START_REF] Fayazbakhsh | Less pain, most of the gain: Incrementally deployable icn[END_REF], [START_REF] Kurose | Information-centric networking: The evolution from circuits to packets to content[END_REF]. Compared to ICNs, the design of content-centric architectures has a more practical purpose as they attempt to meet the dataintensive applications needs by improving content delivery performance and reduce traffic overhead.

The high point of the regain of interest in content-centric networking came with the appearance of the CCN proposal [START_REF] Jacobson | Networking named content[END_REF]. Instead of replacing IP in a clean-state approach, the guiding principles of CCN design follow a more conventional approach that aims to show that removing addresses from IP does not call for revisiting IP's engineering principles.

The key objective behind CCN design is to investigate whether the Internet can evolve with few adjustments from a design paradigm that have been exclusively host-centric to a contentcentric paradigm. As a result, the key choices that have driven the design of CCN were taken to enable the reuse of welltried mechanisms and techniques directly borrowed from IP. Examples of the latters include DNS naming syntax, CIDR prefix aggregation, and longest prefix match forwarding.

In this paper, we focus on the communication functionalities of content-centric networks. More especially, we propose the design of a routing architecture that fully benefits from the widespread availability of content in today's large scale networks. Replication or caching have long become common practice amongst service providers who can improve access efficiency by minimizing retrieval time while reducing bandwidth costs.

The main goal of our work is to complement the existing communication framework of CCN with a scalable interdomain routing protocol able to discover routes that allow users to retrieve replicated content from the most efficient location. We adopt an Interest-based approach for content discovery supported by the caching capabilities of content routers. To handle content replication, we borrow from the literature [START_REF] Walfish | Untangling the web from dns[END_REF] the idea of using flat unstructured identifiers in replacement of CCN DNS-based names to avoid the overloaded semantics of URLs.

In summary, the contributions of our work are as follows:

• Interdomain content routing We describe a scalable discovery scheme of routes towards the most efficient replicas for a content item; • Route caching We show the benefit of using the caching capabilities of content routers to store routes discovered for replicated content; • Flat Namespace We introduce an unstructured locationindependent naming space compliant with the replication or migration of content. The remainder of this paper is structured as follows. In Section II, we present the main features of CCN with regards to naming, routing, and forwarding. Section III provides the detailed description of our routing architecture including the namespace, the node model as well as the types of packets, and finally our routing algorithm. Results of the simulation experiments are presented in Section IV. Finally, the paper is concluded by Section V.

II. CCN ROUTING FRAMEWORK

In this section, we first present the main CCN architectural features with regards to naming, routing, and forwarding. We then draw on the basic principles underlying the design of our routing solution for content centric networks.

1) CCN namespace: CCN name syntax is inspired by the naming scheme of DNS. A CCN name is designed as a URL-like chain of characters with the aim of applying prefix aggregation to the hierarchical structure of domain names. However, the hierarchical structure of DNS names is motivated by the need to distribute the responsibility of assigning domain names. With the current lack of specific allocation policy and management for CCN namespace, the reuse of DNS-based names is expected to result in a mismatch between the name distribution and the underlying network topology which will defy the benefits of hierarchical aggregation. The use of URLlike names also has the drawback of making content replication and migration challenging. In our work, we adopt a flat namespace consisting of unstructured, location-independent, and human-unreadable identifiers to avoid broken links and the use of HTTP redirects which requires the cooperation of the originating administrative domain.

2) Content names distribution: As depicted in Fig. 1, CCN uses the standard OSPF link-state database flooding-based update mechanism to advertise CCN content name prefixes, throughout the routing domain. OSPF has been extended with opaque LSAs that include a field for some application specificinformation. CCN uses the opaque capabilities of OSPF to distribute CCN-specific information i.e. content name prefixes. The CCN routing application attached to the routers can use the opaque LSAs to transparently advertise the content name prefixes throughout the routing domain. The distribution process of CCN-specific information is transparent to the OSPF routing process but can be used by the CCN routing application. Though link-state flooding has been proved to improve routing convergence in medium to large size enterprise or campus networks, OSPF cannot be extended so as to provide interdomain connectivity in large-scale networks. In our work, we leverage the multiprotocol extensions of BGP (i.e. MP-BGP) to carry routing information for the CCN layer.

3) CCN forwarding: CCN applies the concept of aggregation introduced with CIDR addressing to the context of URL-based content names. From the use of aggregation as an already well-tried technique, CCN forwarding is expected to achieve scalability. Introduced in the Internet to efficiently manage the limited addressing space and control the growth in size of the forwarding tables, CIDR address allocation scheme achieves scalability by complying with the location of nodes in the network topology. However, CIDR-like aggregation has the drawback of requiring routers to perform longest-prefixmatch forwarding which makes of fast lookup schemes a challenging design issue for CCN routers. Another issue for CCN forwarding concerns how content items are selected when the requested content is replicated at several locations. Careful attention has to be given to the selection of the paths provided in the FIB to each replica. A suggested strategy consists in broadcasting the request to all replicas which poses the problem of receiving multiple copies of the same content item. A design objective common to most content-centric architectures is to improve content delivery performance and reduce traffic overhead. To achieve this objective, we claim that content-based routing should take benefit from the widespread availability of content in large scale networks as data objects may be replicated and cached at several places. By giving emphasis on the content rather than the location, contentcentric networking renders binding of content to endpoints obsolete when providing access to the content. The paradigm underlying content-centric networking thus gives a natural way to make effective use of content replication since end-to-end connectivity to the source of a replica is not required prior to access the data. As such, content-centric networks naturally provide network-level support for discovery of replicated data items.

We adopt an on-demand interdomain flat name-based routing protocol that integrates naming to routing for effective location of content replicas based on network metrics while providing resilience to network failures or content migration. Our approach overcomes the issue of scalability inherent to routing on flat names by making use of the caching capabilities of content routers for storing previously learnt routes.

The remainder of this section details how our architecture works. We first present our content identifier namespace, the node model as well as the data structures and the types of messages, and then our routing algorithm.

A. Content Identifier CID

We introduce a flat namespace consisting of unstructured, location-independent, and human-unreadable identifiers called Content Identifiers (CIDs). The CIDs are said to be semanticfree since they provide the ability to refer to content and services independently of the hosting endpoint regardless of its administrative domain, network location, or the network topology. In contrast to the CCN DNS-based names, content can be replicated or migrate between hosts while avoiding broken links or using HTTP redirection which requires control of former domain, increases traffic, and latency.

The value of CIDs might be a 128-bit long number chosen randomly or could be derived from the cryptographic hash of the public key of the content owner and the content itself for authentication and integrity respectively. The granularity of naming is flexible since a name can refer to a host, a content within it, or at finer granularity any content item with no reference to the hosting endpoint. At the user-level, it is assumed that a mapping service provides applications with the CIDs corresponding to human readable descriptors such as search keywords which in this case will require the mapping service to be carried out by search engines.

B. Message Types

CCN has two types of messages: Interest and Data messages. The Interest message is used to request data by name. The Data message is used to supply data. In CCN, the Interest message contains the DNS-based name prefix of the content item to be retrieved while the requested content is carried in a Data message identified by a name matching the prefix provided in the Interest. In our routing scheme, the Interest and Data messages now contain a CID field instead of the CCN DNS-based names.

We also extend both types of messages with a field borrowed to BGP update messages. This new field contains a pathvector similar to the BGP AS Path attribute which includes the list of Autonomous System Numbers that describes the sequence of ASes through which the Interest has passed. The path-vector is used to prevent Interests from looping and is returned to the querier in the Data message carrying the requested content item. All subsequent Interests will be directed towards this destination CR of the routing domain where the content was found. The path-vector is also used by content routers to aggregate the Interest messages requesting the same content. A content router stores the shortest path vector followed by previous occurrences of an Interest and only forwards a subsequent one if received along a shorter route. An Interest message is also scoped by a TTL which is initialized to the maximum number of ASes the Interest message can visit.

C. Data Structures

We complete the operation of a CCN node by introducing an additional data structure to the proposed node model. In addition to the three existing structures, namely the FIB (Forwarding Information Base), the Content Store (buffer memory), and the PIT (Pending Interest Table ), we introduce a Route Store (RS). The FIB is used to forward Interests toward the sources having a copy matching the requested content. The Content Store is a cache where content objects can be stored to satisfy future Interests. The PIT keeps track of the interfaces on which Interests for a specific content are forwarded. We extend the PIT entries with information collected from the Interest header fields.

The RS is used to keep in cache the routes discovered in response to previous Interests. Each RS entry contains a pathvector representing the list of ASes that the route traverses.

When a content router (CR) starts searching the best content replica with respect to its relative location, the CR sends an Interest either along the unicast path provided by BGP or will alternatively broadcast the Interest to all its peers in the neighboring ASes depending on the popularity of the content. If a route for the requested content is present in the RS, the CR will use this route to direct the Interest towards the domain where a replica of the content is known to be found.

We extend the PIT entries in order to avoid propagating multiple instances of the same Interest. Each PIT entry is extended with the identifier of the CR who initiated the search, the sequence number of the Interest, and the shortest path vector followed by previous occurrences of this Interest sharing the same originating CR and sequence number. As in CCN, the PIT is also used to aggregate multiple Data messages so as only a single copy will be forwarded towards the content requester.

D. Routing on Content IDs

In this section, we describe how the design of our namebased interdomain routing protocol supports the efficient discovery of some piece of content by using the plurality of alternative routes to access replicated content.

In contrast to CCN, Interests are propagated between neighboring domains in hope of finding the closest content replica. Our choice to avoid designing an advertisement-based protocol is motivated by the need of resiliency to content migration and network failures. Relying on periodic advertisements of routing information is expected to increase the impact on the network when transposed to highly replicated content.

1) Domain: We will refer to a domain as the administrative entity corresponding to a routing domain or an Autonomous System. The boundary of a domain may also be defined by reflecting other organizational considerations related to the ownership of content available within the domain boundary. As such a domain may refer to a campus, a company, or a content provider.

A domain distinguishes content depending whether it is the originator of the content or a purveyor of a copy of the content. A domain containing a local copy of a content considers this content as an internal content. A home domain is a specific case since it contains the original copy of a content item provided by the content owner himself/herself. We assume that an intradomain routing protocol creates the routes needed by local content routers to forward the Interests along the shortest path towards an internal content.

A content item is considered as external to a domain if a copy cannot be found locally. In the absence of a local copy, a domain can still have a route previously learnt for this content. If not, a border content router (BCR) will initiate a search by sending the Interest to its neighboring domains in the same way BGP propagates its update messages. Border content routers are configured to be BGP peers and as so, will establish a BGP session running on top of TCP prior to the exchange of the Interest messages throughout the search process. Inside each domain, a BCR is in charge of propagating the Interests requesting an external content.

2) Requesting an external content: As depicted in Fig. 2, when an Interest for a piece of content replicated in other domains is received, the BCR first checks if the path-vector included in the Interest contains its own domain ID in which case the Interest will be discarded. In the absence of loops, the BCR appends the local domain ID (or alternatively, the AS Number or ASN) to the Interest path-vector and checks the PIT for an entry matching the incoming Interest, i.e. containing the same original requester and sequence number. A PIT-match indicates that a previous Interest has already been sent for the same content. The Interest is discarded if the path-vector included in the PIT entry is shorter than the one carried by the Interest.

Otherwise the PIT entry is updated with the Interest pathvector and arrival interface. If no PIT entry is found, a new PIT entry is created with the identifier of the BCR who initiated the search, the sequence number of the Interest, the path-vector it carries, and the face on which the Interest arrived. The BCR then looks in its Content Store for a matching content. If found, a Data message carrying the content will be returned to the requester by following the PIT entries left behind the Interest.

In case of no Content Store match, the BCR checks in its Route Store (RS) for a route recently learnt through a previous search. If such a route exists, the Interest is forwarded along this route towards the content. If no matching route can be found in the RS, the BCR decrements the TTL of the Interest message. If the Interest TTL did not reach zero, the BCR sends out the Interest to all its neighboring peers or along the unicast BGP route.

A timer is associated to each PIT entry to limit the rate of Interests and to handle unsuccessful searches. If no matching Data is received, the original requester is responsible for generating a new Interest distinguishable from previous occurences thanks to an incremented sequence number.

3) Responding with an internal content: On receiving an Interest for an internal content, a BCR starts performing loopdetection on the Interest path vector in a similar way to BGP. The BCR will append the ID of its domain to the path-vector if loop-free. The resulting sequence of domains corresponds to the complete route between the originating BCR that initiated the search and the local replica of the content. This route is sent back to the originating BCR along with the copy of the local replica in a Data message. The Data message also inherits the Interest sequence number and includes the identifier of the replying BCR.

Upon receiving the Data message, a content router extracts the path vector contained in its header and stores the corresponding route in its Route Store along with the identifier of the BCR who originated the Data message. This path vector will be used to direct subsequent Interests requesting the same content. The CR also stores the content object carried by the Data message in its Content Store before forwarding the Data message on every interfaces listed in the PIT entry left by previous Interests. The PIT entry is erased once the Data The route contained in the Data message will be advertised within the domain of the original requester so that internal content routers may also learn how to reach the external content in case of subsequent Interests. This route will also be used so as to handle Interests received from external peer domains.

4) Default routing for sparse content: Since some content can exhibit a low degree of replication, we introduce an alternative using the unicast route towards the home domain of the requested content. In order to enable the use of CIDs for unicasting Interests, we use the multiprotocol extensions of BGP (i.e. MP-BGP) [START_REF] Bates | Multiprotocol extensions for bgp-4[END_REF] to carry routing information for the CCN layer using CIDs.

We propose to define a new address family called CCN AF used to identify CCN as the network layer protocol carried by BGP. To make BGP advertise a set of CIDs as reachable destinations together with the next hop identifier to be used for forwarding along the feasible routes to these destinations, we use the CCN AF in combinaison with the two BGP attributes Multiprotocol Reachable NLRI (MP_REACH_NLRI) and Multiprotocol Unreachable NLRI In our protocol, the Address Family information identifies CCN as the network layer protocol that is being carried within the update. The Next-Hop Information is the identifier of the next content border router on the path to the destination. In our protocol, the NLRI is encoded as one or more 2-tuples with the following format: Length (1 octet): Total length of the CID advertised; CID: The content identifier.

IV. EVALUATION

A. Experimental Setup

In order to evaluate our proposal, we developed ccnSim++, an extension to the ccnSim simulator [START_REF] Chiocchetti | ccnsim: An highly scalable ccn simulator[END_REF]. ccnSim is a chunklevel simulator built in C++ and developed on top of the Omnet++. Its ability to scale well in dense topologies allowed us to conduct our simulations on a large dataset made publicly available by the Cooperative Association for Internet Data Analysis (CAIDA). This dataset consists of a snapshot of more than 36.000 ASes and around six million announced BGP routes. In our simulations, we assume that all links have unlimited bandwidth and are not affected by congestion. All our tests were performed on a machine equipped with 144GB of RAM and a 24 core Intel Xeon processor.

Inside each Autonomous System, we suppose all content objects are forwarded efficiently according to an intradomain routing protocol. Border routers are considered to be content repositories and are responsible for handling Interest and Data messages. To simulate content placement and replication over the Internet, we distribute content replicas randomly across all repositories. The catalog we consider consists of 10 4 unique content objects.

ccnSim++ generates Interests messages according to a Zipf distribution. By keeping track of the number of Interests, every router is able to decide locally which content is popular. As better routes to closer content are only searched for popular content, we simulate the routes lifetime by implementing a LRU Route Store. When the cache fills up, the routes to repositories containing less popular content are removed and replaced with new ones. Fig. 3 shows the length ratio of the BGP route to the original content over the routes discovered by our solution to a replica as a function of the replication degree of the content considered. We can see that the routes we discover outperform the routes to the original content as provided by BGP. This result indicates that our discovery scheme is able to always find the shortest path to the closest content replica. BGP has no way of detecting the repositories holding a replica of the requested content. Moreover, BGP policies may prevent the selection of the shortest path among the routes advertised toward the original content. The length ratio converges to 1 as the number of replicas increases. Due to the random distribution of all content replicas, the placement of replicas is not guaranteed as optimal with regard to the queriers. For large number of replicas, the distant replicas outbalance the replicas close to each querier. The BGP route provided for the original content shares similarities with the routes discovered to the replicas as their number increases and more replicas are well scattered across the network. In all other cases, our query-based content discovery scheme outperforms the BGP path.

C. The Influence of Popularity

In figure Fig. 4, we show the percentage of routes that are actually discovered when searching for a content or reused from a previous successful search as a function of the popularity of the content which follows a Zipf distribution. We can see that over than 75% of the routes used to retrieve a content are provided from the Route Store of the content routers. This shows the benefit of caching the routes discovered toward a replica as a result of requesting a content. The overhead of the discovery process is outweighed by the benefit of caching the routes to the replicas we discover. The benefit of caching routes toward previously discovered replicas can also be found in the number of Interest messages. Fig. 5 shows the total number of Interest messages forwarded in the network as a function of the popularity of the content. Though the number of Interest messages issued for a content increases with the popularity of the content, we can see an opposite behavior when we increase the Zipf skew parameter. This indicates that Interest messages travel less distance in the network to discover a content as its popularity increases. Interest messages can draw on previous searches as content routers intercept them and reply with cached routes With a large catalog size, the route store is unable to hold all the best routes to the popular content. In Fig. 6, we evaluate the impact of the size of the content catalogue as we expect the limited storage of Route Stores to impact of the performance of our route discovery scheme. We vary the size of the Route Stores for different sizes of the content catalogue. We measure the resulting number of Interest and Data messages as a function of the ratio of the RS size over the catalogue. We can observe that for small increases in the size of the Route Stores, the number of messages generated by our architecture decreases significantly. For size ratios greater than 0.8%, this number of Interest and Data messages remains stable. This is due to the fact that the size of the route stores becomes large enough in comparison to the catalogue size to hold the routes discovered for all popular contents.

V. CONCLUSION

In this paper, we have presented a scalable on-demand routing protocol that takes advantage of the paradigm of content-centric networking to meet the needs of content related services and data-intensive applications. As part of our design process, we take advantage of the high availability of content in large-scale networks. The novelty of our approach lies in the use of the caching capabilities of content routers to store routes discovered for content replicas. Our protocol complements the CCN routing framework with the ability to discover efficient routes that allow users to retrieve replicas in a scalable fashion.

Our simulation results show an increase in access efficiency by minimizing retrieval time and by maximizing the bandwidth demand. As a future work, we aim to study the resilience of our protocol to the episodic availability of hosting nodes due to network failures or content migration.

Fig. 1 :

 1 Fig. 1: CCN communication framework

Fig. 2 :

 2 Fig. 2: Handling Interests

Fig. 3 :

 3 Fig. 3: Discovered route evaluation

Fig. 4 :

 4 Fig. 4: Discovery process benefit

Fig. 5 :Fig. 6 :

 56 Fig. 5: Search messages overhead

ACKNOWLEDGMENT

This work was partially funded by the French National Research Agency (ANR) under project ANR DataTweet (ANR-13-INFR-0008).