Soumission

Thomas Opitz
email: thomas.opitz@paca.inra.fr

Latent Gaussian modeling and INLA: A review with focus on space-time applications Titre: Modèles à processus gaussiens latents et inférence INLA : un survol orienté vers les applications spatio-temporelles

Keywords: Integrated Nested Laplace Approximation, R-INLA, spatio-temporal statistics Mots-clés : Integrated Nested Laplace Approximation, R-INLA, statistique spatio-temporelle

Bayesian hierarchical models with latent Gaussian layers have proven very flexible in capturing complex stochastic behavior and hierarchical structures in high-dimensional spatial and spatio-temporal data. Whereas simulation-based Bayesian inference through Markov Chain Monte Carlo may be hampered by slow convergence and numerical instabilities, the inferential framework of Integrated Nested Laplace Approximation (INLA) is capable to provide accurate and relatively fast analytical approximations to posterior quantities of interest. It heavily relies on the use of Gauss-Markov dependence structures to avoid the numerical bottleneck of high-dimensional nonsparse matrix computations. With a view towards space-time applications, we here review the principal theoretical concepts, model classes and inference tools within the INLA framework. Important building blocks for space-time models are certain spatial Matérn-like Gauss-Markov random fields, obtained as approximate solutions to a stochastic partial differential equation. Efficient implementation of statistical inference tools for a large variety of models is available through the INLA package of the R software. To showcase the practical use of R-INLA and to illustrate its principal commands and syntax, a comprehensive simulation experiment is presented using simulated non Gaussian space-time count data with a first-order autoregressive dependence structure in time.

Résumé : Les modèles bayésiens hiérarchiques structurés par un processus gaussien latent sont largement utilisés dans la pratique statistique pour caractériser des comportements stochastiques complexes et des structures hiérarchiques dans les données en grande dimension, souvent spatiales ou spatio-temporelles. Si des méthodes d'inférence bayésienne de type MCMC, basées sur la simulation de la loi a posteriori, sont souvent entravées par une covergeance lente et des instabilités numériques, l'approche inférentielle par INLA ("Integrated Nested Laplace Approximation") utilise des approximations analytiques, souvent très précises et relativement rapides, afin de calculer des quantités liées aux lois a posteriori d'intérêt. Cette technique s'appuie fortement sur des structures de dépendance de type Gauss-Markov afin d'éviter des difficultés numériques dans les computations matricielles en grande dimension. En mettant l'accent sur les applications spatio-temporelles, nous discutons ici les principales notions théoriques, les classes de modèles accessibles et les outils d'inférence dans le contexte d'INLA. Certains champs Gauss-Markov spatiaux, obtenues comme solution approximative d'une équation différentielle partielle stochastique, sont des briques de construction importantes dans la modélisation spatio-temporelle. Pour illustrer l'utilisation pratique du logiciel R-INLA et la syntaxe de ses commandes principales, un scénario de simulation-réestimation est présenté en détail, basé sur des données simulées, spatio-temporelles et non gaussiennes, avec une structure de dépendance autorégressive dans le temps.

Introduction

The rapidly increasing availability of massive sets of georeferenced data has spawned a strong demand for suitable statistical modeling approaches to handle large and complex data. Bayesian hierarchical models have become a key tool for capturing and explaining complex stochastic structures in spatial or spatio-temporal processes. Many of these models are based on latent Gaussian processes, embedded in a parameter related to the mean of the distribution assumed for the likelihood of the data, and extend the Gaussian random field modeling brought forward by classical geostatistics. Using a conditional independence assumption for the data process with respect to the latent Gaussian layer makes inference tractable in many cases. Typically, closed-form expressions for the likelihood are not available for these complex models, and simulation-based inference through Markov chain Monte Carlo has become a standard approach for many models. An important alternative, superior to MCMC inference under certain aspects, has been developed through the idea of Integrated Nested Laplace Approximation, proposed in the JRSS discussion paper of [START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF]. Many case studies have been conducted through INLA in the meantime, with space-time applications to global climate data [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF], epidemiology [START_REF] Bisanzio | Spatio-temporal patterns of distribution of West Nile virus vectors in eastern Piedmont Region, Italy[END_REF], disease mapping and spread [START_REF] Schrödle | Spatio-temporal disease mapping using INLA[END_REF][START_REF] Schrödle | Assessing the impact of a movement network on the spatiotemporal spread of infectious diseases[END_REF], forest fires [START_REF] Serra | Spatio-temporal log-Gaussian Cox processes for modelling wildfire occurrence: the case of Catalonia, 1994-2008[END_REF], air pollution risk mapping [START_REF] Cameletti | Spatio-temporal modeling of particulate matter concentration through the SPDE approach[END_REF] or fishing practices [START_REF] Cosandey-Godin | Applying Bayesian spatiotemporal models to fisheries bycatch in the Canadian Arctic[END_REF]. More generally, INLA has been successfully applied to generalized linear mixed models [START_REF] Fong | Bayesian inference for generalized linear mixed models[END_REF], log-Gaussian Cox processes [START_REF] Illian | A toolbox for fitting complex spatial point process models using integrated nested Laplace approximation (INLA)[END_REF], survival models [START_REF] Martino | Approximate Bayesian inference for survival models[END_REF] and quantile regression [START_REF] Yue | Bayesian inference for additive mixed quantile regression models[END_REF], amongst many other application fields. The recent monograph of [START_REF] Blangiardo | Spatial and Spatio-temporal Bayesian Models with R-INLA[END_REF] reviews INLA in detail and gives many practical examples. Instead of applying simulation techniques to produce a representative sample of the posterior distribution, INLA uses analytic Laplace approximation and efficient numerical integration schemes to achieve highly accurate analytical approximation of posterior quantities of interest with relatively small computing times. In particular, we get approximations of univariate posterior marginals of model hyperparameters and of the latent Gaussian variables. By making use of latent Gauss-Markov dependence structures, models remain tractable even in scenarios that are very high-dimensional in terms of observed data and latent Gaussian variables.

The INLA-based inference procedures are implemented in the R-package INLA (referred to as R-INLA in the following) for a large variety of models, defined through basic building blocks of three categories: the (univariate) likelihood specification of data, the latent Gaussian model and prior distributions for hyperparameters. Functionality of R-INLA is continuously extended [START_REF] Martins | Bayesian computing with INLA: new features[END_REF][START_REF] Lindgren | Bayesian spatial modelling with R-INLA[END_REF][START_REF] Rue | Bayesian computing with INLA: A review[END_REF]. This review and the code examples refer to R-INLA version 0.0-1463562937. The R-INLA software project is hosted on http://www.r-inla.org/, where one can find lots of INLA-related resources, amongst them details on the specification of likelihoods, latent models and priors, a discussion forum with very active participation of the members of the INLA core team, tutorials and codes, an FAQ section, etc.

Modeling and estimation with INLA

Latent Gaussian modeling

The structured latent Gaussian regression models amenable to INLA-based inference can be defined in terms of three layers: hyperparameters, latent Gaussian field, likelihood model. The univariate likelihood captures the marginal distribution of data and is usually chosen as an exponential family (Gaussian, gamma, exponential, Weibull, Cox, binomial, Poisson, negative binomial, ...), where models like the exponential, Weibull or Cox ones are available as survival models allowing for right-and left-side censoring. Following the mechanism of generalized linear models, the mean of the likelihood distribution is determined by the latent Gaussian predictor through a link function such that E(Y | η) = h -1 (η), where Y is the observation, η is a Gaussian predictor and h is an appropriately chosen link function. Hyperparameters can appear in the likelihood as dispersion parameters like the variance of the Gaussian distribution, the overdispersion parameter of the negative binomial one or the shape parameter of the gamma one, or they can characterize the structure of the latent Gaussian model, for instance through variances, spatial correlation parameters or autoregression coefficients. Formally, this hierarchical model can be written as θ θ θ ∼ π(θ θ θ ) hyperparameters (1)

x x x | θ θ θ ∼ N (0 0 0, Q Q Q(θ θ θ ) -1 ) latent Gaussian field (2) y y y | x x x, θ θ θ ∼ ∏ i π(y i | η i (x x x), θ θ θ ) observations (3)
where Q Q Q(θ θ θ ) is the precision matrix (i.e., inverse covariance matrix) of the latent Gaussian vector x x x, η η η(x x x) = A A Ax x x with the so-called observation matrix A A A that maps the latent variable vector x x x to the predictors η i = η i (x x x) associated to observations y i . If y y y and x x x can be high-dimensional when using INLA, an important limitation concerns the hyperparameter vector θ θ θ ∈ Θ whose dimension should be moderate, say ≤ 10, which is due to numerical integration that has to be carried out over the hyperparameter space Θ. Notice that the Gaussian likelihood is particular since, conditional to the hyperparameters, the observations are still Gaussian. In practice, the precision hyperparameter of the Gaussian likelihood (i.e., the inverse of its variance) can correspond to a measurement error or a nugget effect, and we can fix a very high value for the precision hyperparameter if we want the model for data y y y to correspond exactly to the latent Gaussian predictors η η η. The dependence structure between observations y i is captured principally by the precision matrix Q Q Q(θ θ θ ) of the latent field x x x. In practice, it is strongly recommended or even indispensable from the point of view of computation time and memory requirements to choose Gauss-Markov structures with sparse Q Q Q(θ θ θ ) whenever model dimension is high.

The resulting joint posterior density of latent variables x x x and hyperparameters θ θ θ is

π(x x x, θ θ θ | y y y) ∝ exp -0.5x x x Q Q Q(θ θ θ )x x x + ∑ i log π(y i | η i , θ θ θ ) + log π(θ θ θ ) . (4)
This density over a high-dimensional space does usually not characterize one of the standard multivariate families and is therefore difficult to interpret and to manipulate. In practice, the main interest lies in the marginal posteriors of hyperparameters θ j , of latent variables x i and of the resulting predictors η i , which can be included into x x x for notational convenience. Calculation of these univariate posterior densities requires integrating with respect to θ θ θ and x x x:

π(θ j | y y y) = π(x x x, θ θ θ | y y y)dx x x dθ θ θ -j = π(θ θ θ | y y y)dθ θ θ -j , (5)
π(x i | y y y) = π(x x x, θ θ θ | y y y)dx x x -i dθ θ θ = π(x i | θ θ θ , y y y)π(θ θ θ | y y y) dθ θ θ . ( 6
)
We notice that the use of astutely designed numerical integration schemes with respect to the moderately dimensioned hyperparameter space Θ can yield satisfactorily accurate approximation of the outer integral. On the other hand, calculating the inner integral with respect to x x x, often of very high dimension (≈ 10 2 to 10 8 ), is intricate.

Gauss-Markov models

We say that a random vector

x x x | θ θ θ ∼ N (0 0 0, Q Q Q -1 ) is Gauss-Markov if the number of nonnull entries of its n × n precision matrix Q Q Q is O(n).
Such sparse precision matrices allow efficient numerical computation of matrix operations like LR-decomposition (with sparse factors L and R), determinant calculation, matrix-vector products, etc. For instance, complexity of matrix inversion decreases from O(n 3 ) for matrices without any structural constraints to around O(n 3/2 ) for sparse matrices. Using Gauss-Markov structures fundamentally shifts the dependence characterization from covariance matrices

Q Q Q -1 to precision matrices Q Q Q. Notice that the conditional expectation is easily expressed through the regression E(x i | x x x -i ) = -∑ j =i (q i j /q ii )x j
, where only a small number of the sum terms, also called the neighborhood of x i , are non-zero owing to the sparse structure of

Q Q Q. The conditional variance is V(x i | x x x -i ) = 1/q ii . R-INLA
uses fast and efficient algorithms for sparse matrix calculations [START_REF] Rue | Gaussian Markov random fields: theory and applications[END_REF], already implemented in the GMRFLib library. For efficient calculations, it is important to make the precision matrix Q Q Q "as diagonal as possible" by reordering variables to regroup nonzero elements as close as possible to the diagonal. R-INLA has implemented several of those reordering strategies; see [START_REF] Rue | Gaussian Markov random fields: theory and applications[END_REF] for more details on reordering algorithms. If certain Gauss-Markov models exist for spatially indexed graphs, useful covariance functions defined over R d and leading to Gauss-Markov covariance matrices are difficult to establish. An exception is the very flexible approximate Gauss-Markov representation of Matérn-like covariances based on certain stochastic partial differential equations (often referred to as the SPDE approach in the literature), which is also implemented in R-INLA; see Section 3.1 for more details.

INLA

The fundamental idea of INLA consists in applying the device of Laplace approximation to integrate out high-dimensional latent components. This theoretical foundation is combined with efficient algorithms and numerical tricks and approximations to ensure a fast yet accurate approximation of posterior marginal densities of interest like those of the predictors η i in (6) or of hyperparameters θ j in (5). Since the details of methods implemented in the INLA approximation are quite technical, we here content ourselves with a presentation of the main ideas and the related options available in R-INLA.

The principle of Laplace approximation

We first recall the principle of the Laplace approximation and its calculation in practice. Typically, one seeks to evaluate an integral f (x x x) dx x x, where the positive integrand function f , here written as f (x x x) = exp(g(x x x)), is defined over a high-dimensional space and is "well-behaved" in the sense that it is unimodal and that its shape is not too far from gaussianity. Since the integral value is mainly determined by the behavior around the mode of g, a second-order Taylor approximation of g can be substituted for g to calculate an approximate value of the integral. Assuming that x x x is the unique global maximum of g, we get g(x x x) ≈ g(x x x ) + 0.5(x x xx x x ) H H H(g)(x x x )(x x xx x x ) for values x x x close to x x x with the Hessian matrix H H H(g)(x x x ). Notice that -H H H(g)(x x x ) is positive definite. More generally, we now write f (x x x) = exp(kg(x x x)) with a scale variable k > 0. An approximate value of the integral can be calculated using the fact that a multivariate Gaussian density integrates to 1.

The resulting following integral approximation in dimension d is expected to become more and more accurate for higher values of k 1:

∞ ∞ ∞ -∞ ∞ ∞ f (x x x) dx x x = ∞ ∞ ∞ -∞ ∞ ∞ exp(kg(x x x)) dx x x (7) k 0 ≈ ∞ ∞ ∞ -∞ ∞ ∞ exp(kg(x x x ) + 0.5k(x x x -x x x ) H H H(g)(x x x )(x x x -x x x )) dx x x = 2π k d/2 |H H H(g)(x x x )| -1/2 exp(kg(x x x )). (8)
In statistical practice, k may represent the number of i.i.d. replications, each of which has density f (x x x) = exp(g(x x x)). Higher values of k usually lead to better approximation, and formal results on the quality of approximation have been derived [START_REF] Tierney | Accurate approximations for posterior moments and marginal densities[END_REF][START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF]. Many of the models commonly estimated with INLA have no structure of i.i.d. replication, but the Laplace approximation remains sufficiently accurate in most cases since there usually still is a structure of internal replication, for instance when there are blocks of data that are close to independence.

In the context of INLA, the following observation will be interesting and useful. Fix k = 1 in (7) and suppose that f (x x x) = exp(g(x x x)) = π(x x x, θ θ θ ), where π(x x x, θ θ θ ) is the joint probability density of a random vector (x x x, θ θ θ ). Then, in (8), the term exp(g(x x x )) is the value of π at its mode x x x for fixed θ θ θ , whereas (2π

) d/2 |H H H(g)(x x x )| -1/2 is 1/π G (x x x | θ θ θ
) with π G a Gaussian approximation with mean vector x x x to the conditional density of x x x | θ θ θ . In practice, we can determine the mean µ µ µ = x x x and the precision matrix Q Q Q = -H H H(g)(x x x ) of π G through an iterative Newton-Raphson optimisation. Starting from the joint posterior (4) of our latent Gaussian model, we set g

(x x x) = -0.5x x x Q Q Q(θ θ θ )x x x + ∑ i log π(y i | η i , θ ). We further write g i (x i ) = log π(y i | x i , θ θ θ ) and calculate its second-order Taylor expansion g i (x i ) ≈ g i (µ (0) i ) + b i x i -0.5c i x 2
i . Without loss of generality, we here assume that the linear predictor η η η corresponds to the latent Gaussian vector x x x. We start the iterative optimisation with initial values

Q Q Q (1) = Q Q Q + diag(c c c) and µ µ µ (1) , where Q Q Q (1) µ µ µ (1) = b b b. We then iterate this procedure until convergence such that µ µ µ (i) → µ µ µ = x x x and Q Q Q (i) → Q Q Q = Q Q Q + diag(c c c )
, where an appropriate convergence criterion must be used. Notice that the conditional independence assumption of observations y i with respect to (η i , θ θ θ ) allows preserving the sparse structure in Q Q Q . It is further possible to impose linear constraints M M Mx x x = e e e Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016 onto x x x and x x x with given matrix M M M and vector e e e by using the approach of conditioning through kriging [START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF].

Posterior marginal densities of hyperparameters

To calculate

π(θ j | y y y) = π(x x x, θ θ θ | y y y)dx x x dθ θ θ -j = π(θ θ θ | y y y)dθ θ θ -j , (9)
we use the Laplace approximation of the inner integral π(x x x, θ θ θ | y y y)dx x x = π(θ θ θ | y y y) as described in Section 2.3.1 such that the approximated density π satisfies

π(θ θ θ | y y y) ∝ π(x x x, θ θ θ , y y y) π G (x x x | θ θ θ , y y y) | x x x=x x x (θ θ θ ) (10)
with x x x (θ θ θ ) the mode of the joint density π(x x x, θ θ θ , y y y) for fixed (θ θ θ , y y y) and a Gaussian density π G that approximates π(x x x | θ θ θ , y y y):

π G (x x x | θ θ θ , y y y) = (2π) n/2 |Q Q Q (θ θ θ )| 1/2 exp -0.5(x x x -x x x (θ θ θ )) Q Q Q (θ θ θ )(x x x -x x x (θ θ θ )) . (11)
Notice that the Gaussian approximation π G is exact if the data likelihood π(

y i | η i , θ θ θ ) itself is
Gaussian. An approximation of the posterior marginal of θ j in ( 9) is now obtained through a numerical integration with a set of integration nodes θ θ θ chosen from a numerical exploration of the surface of the density π(θ θ θ -j , θ j | y y y) (with θ j held fixed). This yields

π(θ j | y y y) = L ∑ =1 ω π(θ θ θ | y y y) (12)
with weights ω (which are chosen to be 1 in the approaches implemented in R-INLA). In R-INLA, θ θ θ can either be chosen as a grid around the mode of π(θ θ θ | y y y) (int.strategy="grid", the most costly variant), or through a simpler so-called complete composite design which is less costly when the dimension of θ θ θ is relatively large (int.strategy="ccd", the default approach), or we may use only one integration node given as the mode value (int.strategy="eb", corresponding to the idea of an empirical Bayes approach).

Posterior marginal densities of the latent Gaussian field

For calculating the marginal density π(x i | y y y) of a latent variable x i , we lean on representation (6). Numerical integration with respect to θ θ θ can be done in analogy to the procedure described in Section 2.3.2, and the Laplace approximation (10) allows approximating π(θ θ θ | y y y). It thus remains to (approximately) evaluate π(x i | θ θ θ , y y y). A simple and fast solution would be to use the univariate Gaussian approximation resulting from the multivariate Gaussian approximation ( 11) whose mean value is x i (θ θ θ ) and whose variance can easily and quickly be calculated from a partial inversion of the precision

Q Q Q (θ θ θ ) (Rue, 2005) (strategy="gaussian" in R-INLA).
However, this Gaussian approximation often fails to capture skewness behavior and can generate nonnegligible bias in certain cases -an important exception to this issue being the case where the data likelihood is Gaussian. In the general case, using again a Laplace-like approximation

π(x x x, θ θ θ , y y y) π G (x x x -i | x i , θ θ θ , y y y) | x x x -i =x x x -i (x i ,θ θ θ ) (13)
with mode x x x -i (x i , θ θ θ ) of π(x x x, θ θ θ , y y y) for fixed (x i , θ θ θ , y y y) would be preferable, but is relatively costly (strategy="laplace" in R-INLA). Instead, [START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF] propose a so-called simplified Laplace approximation based on third-order Taylor developments of numerator and denominator in ( 13) that satisfactorily remedies location and skewness inaccuracies of the Gaussian approximation (strategy="simplified.laplace" in R-INLA, the default). Notice that the "Nested" in INLA refers to this second Laplace-like approximation.

Space-time modeling approaches

Modeling nonstationarities over time and spatio-temporal dependence in repeated spatial observations is paramount to understanding the dynamics of processes observed over space and time. We here review approaches to integrating the time component into the latent Gaussian predictor η st that are suitable for high-dimensional space-time inference with INLA. In principle, any space-time Gaussian process η st could be used, but the requirement of a Gauss-Markov structure for fast matrix calculations and the current scope of models implemented in R-INLA impose some constraints. Flexible, Matérn-like spatial Gauss-Markov models with continuous sample paths are available in R-INLA and will be discussed in Section 3.1. A generalization of such purely spatial models to flexible nonseparable space-time dependence structures is still pending, but R-INLA allows extending such spatial models to capture temporal dependence through autoregressive structures, which will be discussed in the following and explored in the examples of Section 4. New approaches to generic nonseparable space-time Gauss-Markov model classes and their implementation in R-INLA are projected by members from the scientific community working on the SPDE approach. Any other nonseparable covariance model would in principle be amenable to INLA-based inference without difficulty as long as it keeps a Markovian structure when using a high-dimensional latent model. We start by formulating a generic latent Gaussian space-time model for the predictor η st that covers many of the models that can be fitted with R-INLA. We denote by ζ stk , k = 1, . . . , K given covariate data, which may depend on space or time only. Notice that in cases where covariates are not available except at the observation points (s i ,t i ), i = 1, . . . , n of y i , one may have to include them into into the observation vector y y y and formulate a latent Gaussian model for their interpolation over space and time. This would give a model that achieves both interpolation of covariate values and prediction of y y y at the same time. Here, we consider

η st = x 0 + K ∑ k=1 x k ζ stk + L ∑ k=1 f k (ζ stk ) + x t + x s + x st (14)
The linear coefficients x 0 , . . . , x K are known as fixed effects, whereas functions f k (•) and the processes x t , x s and x st are referred to as random effects. Notice that x t and x s are to be understood as marginal models that capture purely spatial or temporal marginal effects. We now shortly present some typical examples of Gaussian prior models that are commonly used and available in R-INLA for the intercept x 0 , the linear covariate effects x k , the nonlinear covariate effects f k (•), the marginal temporal effect x t , the marginal spatial effect x s and the space-time effect x st .

A temporal effect x t could be modeled as an autoregressive process or as a random walk, where autoregression coefficients and the step variance of the random walk represent hyperparameters that could be estimated. A nonparametric spatial effect x s could be modeled with a Gauss-Markov Matérn-like prior random field, the details of whose construction are presented in the following section 3.1. A simple extension to a space-time effect x st is obtained from considering independent replicates of the spatial field for each time point. An important class of space-time models x st that allows for temporal dependence and preserves the Gauss-Markov structure are the stationary first-order autoregressive models

x st = ax s,t-1 + 1 -a 2 ε st , (15)
where a ∈ (-1, 1] and ε st is a stationary spatial innovation process, i.i.d. in time, typically chosen as the Gauss-Markov Matérn-like field. We here allow for a = 0 to include the case of temporal independence as mentioned before. This AR model is a group model where spatial groups ε st , t = 1, 2, . . ., are joined through an AR group model. In R-INLA, it is possible to work with more general group models that define a type of dependence like "autoregressive", "exchangeable", "random walk" or "besag" between certain groups of latent variables. It is further possible to specify certain graph structures among latent variables; we here refer to www.r-inla.org for full details about the specification of a large variety of available latent models. Owing to issues of identifiability and model complexity, usual only a subset of the terms in ( 14) is used to construct the latent field in practical applications.

Spatial Gauss-Markov random fields based on the SPDE approach

The spatial SPDE model of [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] defines a Gauss-Markov random field as the approximate solution to a certain stochastic partial differential equation. It is an important building block for latent Gaussian models with spatial and spatio-temporal effects. Contrary to classical covariance function models, this approach provides sparse precision matrices that make numerical procedures efficient even for very high-dimensional problems. Formally, a Gaussian process x(s) on R D is defined through

κ 2 -∆ α/2 x(s) = W (s), α = ν + D/2, s ∈ Ω (16)
with the Laplace operator ∆y = ∑ D j=1 ∂ 2 y/∂ 2 x j , a standard Gaussian white noise process W (s) and a nonempty spatial domain Ω ⊂ R D with regular boundary. Depending on the value of ν and D, the Laplace operator κ 2 -∆ α/2 is fractionary with noninteger exponent α/2, and it must be defined in an appropriate way [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF]. The only stationary solution to (16) for Ω = R D is a Gaussian random field with the Matérn covariance function whose shape parameter is ν (with ν = 0.5 yielding the exponential covariance model). The marginal variance is 0.1 is attained between two points is around √ 8ν/κ 2 . The Matérn model is known to be very flexible through its scale and shape parametrization, with regularity properties of sample paths governed by the shape parameter ν.

Γ(ν)/[Γ(ν + D/2)(4π) D/2 κ 2ν ],
In practice, when working on a finite domain Ω ⊂ R 2 , boundary effects come into play. One can assume a polygon-shaped boundary ∂ Ω, as it is implemented in R-INLA. An interesting choice of boundary condition is the Neumann condition with zero normal derivatives at the boundary such that the Gaussian field is "reflected" at the boundary. [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] show that Neumann conditions principally lead to an increase in variance close to the boundary, the factor being approximately 2 when there is one close linear boundary segment, and 4 when we are close to the 90-degree angle of a rectangle where two linear segments meet. Whereas such boundary conditions may be interesting for some applications, we often prefer to extend the domain Ω beyond the study region towards a larger domain, such that boundary effects become negligible within the study region. This requires that the extended domain's boundary is separated by a distance superior to the empirical range from the study region.

Approximate solutions to ( 16) are obtained based on the finite element approach commonly used in the numerical solution of partial differential equations. Using a triangulation of the spatial domain leads to a high-dimensional multivariate representation with Gaussian variables located in the triangulation nodes s i . Spatial interpolation between nodes is achieved by considering these Gaussian variables as weights for piecewise linear basis functions ψ i (s), one for each node. Basis functions are of compact support giving by the triangles touching the node s i ("pyramid functions"); see Figure 1 for an example of a basis function and of the approximation of a spatial surface through a linear combination of such basis functions. By "projecting" the SPDE ( 16) on the space spanned by the Gaussian variables, one can calculate the precision matrix Q Q Q governing the dependence structure between these variables. There are certain rules of thumb to be respected for a construction of the triangulation that does not strongly distort the actual dependence structure of the exact solution to (16) and that remains numerically stable with respect to certain matrix computations, mainly concerning maximum sizes of triangles and minimum sizes of interior angles. For numerical efficiency, overly fine triangulations can be avoided by requiring minimum edge lengths, for instance; see [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] for further details on the finite element construction. The approximate Gauss-Markov solution has mean zero and its precision matrix Q Q Q has entries that are determined by the SPDE. We refer to the Appendix section of [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] for the calculation of Q Q Q, which is explicit. Based on the approximate solutions for α = 0 and α = 1, an approximate solution of the SPDE for α + 2 can be obtained by injecting the solution for α at the place of the white noise W (s) in ( 16). For non-integer values of α > 0, additional approximations are necessary to obtain a solution, see the authors' response in the discussion of the [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] paper. We have here considered κ and τ to be constant over space. It is possible to allow spatial variation of these parameters for nonstationary models with κ = κ(s) as in ( 16) and a precision-related parameter τ = τ(s) that varies over space. To wit, [START_REF] Ingebrigtsen | Spatial models with explanatory variables in the dependence structure[END_REF] apply such second-order nonstationary modeling to precipitation in Norway, where altitude is used as a covariate that acts on the local covariance range in the dependence structure. The SPDE ( 16) is also well-defined over manifolds Ω, for instance the sphere S 2 embedded in R 3 . In general, [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] show that the approximate solution converges to the true solution of the SPDE for an adequately chosen norm when the triangulation is refined in a way such that the maximum diameter of a circle inscribed into one of the triangles tends to 0.

R-INLA currently implements calculation, estimation and simulation of the Gauss-Markov SPDE solution for α ≤ 2, for Ω a subset of R 1 , R 2 or a two-dimensional manifold embedded into R 3 , and for spatially constant or varying values κ(s) and τ(s). A large number of tools is available to construct numerically efficient and stable triangulations.

Using R-INLA on a complex simulated space-time data set

To execute the following code, the INLA package of R must be installed; see www.r-inla.org for information on the R command recommended for installing this package, which is not hosted on the site of the Comprehensive R Archive Network (CRAN) due to its use of external libraries. We here illustrate the powerful estimation and inference tools of R-INLA in a controlled simulation experiment with data simulated from a latent Gaussian space-time model with the Poisson likelihood. The full code for the simulation study below can be requested from the author.

Simulating the data

We simulate a space-time count model based on a latent first-order autoregressive Gaussian process defined on [0, 1] 2 for t = 1, . . . , 60. We use two covariates, given as z 1 (t) = t/60 and z 2 (t) simulated according to an autoregressive time series model. The simulated model is

Y (s,t) | η(s,t) ∼ Pois(exp(η(s,t))) i.i.d. η(s,t) = -1 + z 1 (t) + 0.5z 2 (t) +W (s,t) W (s, 1) = ε(s, 1) W (s,t) = 0.5W (s,t -1) + (1 -0.5) 2 ε(s,t), t = 2, . . . , 60.
The covariance function of the standard Gaussian field ε(s,t) is chosen of Matérn type with shape parameter ν = 1 and effective range 0.25. We now fix 50 × 60 = 3000 observation points (s i ,t i ) of Y (s,t), determined as the Cartesian product of 50 sites uniformly scattered in [0, 1] and t i = i, i = 1, . . . , 60.

To illustrate the simulation capacities of R-INLA, we here use the SPDE approach to achieve simulation based on the Gauss-Markov approximation of the Matérn correlation structure. After loading the INLA-package and fixing a random seed for better reproducibility of results, we start by defining the κ, τ and α parameters of the SPDE. To avoid boundary effects in the SPDE simulation, we wil use the square [-0.5, 1.5] 2 as spatial support. Next, we create a fine 2D triangulation mesh for relatively accurate simulation, with maximum edge length 0.04 within [0, 1] 2 and 0.2 in [-0.5, 1.5] 2 \ [0, 1] 2 . The minimum angle between two edges is set to 21, a value recommended to avoid ill-conditioned triangulations containing very elongated triangles. Polygon nodes should be given in counterclockwise order, whereas as polygon-shaped holes in the support would be specified by clockwise order of nodes; see the left display of Figure 2 for the resulting triangulation.

nodes.bnd=matrix(c(0,0,1,0,1,1,0,1),ncol=2,byrow=T) segm.bnd=inla.mesh.segment(nodes.bnd) nodes.bnd.ext=matrix(c(-.5,-.5,1.5,-.5,1.5,1.5,-.5,1.5),ncol=2,byrow=T) segm.bnd.ext=inla.mesh.segment(nodes.bnd.ext) mesh.sim=inla.mesh.2d(boundary=list(segm.bnd,segm.bnd.ext), max.edge=c(.04,.2),min.angle=21) plot(mesh.sim)

The list slot mesh.sim$n informs us that there are 2401 triangulation nodes. The mesh object mesh.sim has a slot mesh$loc which contains a three column matrix. The first two columns indicate the 2D coordinates of mesh nodes. In our case, the third coordinate, useful for specifying 2D manifolds in 3D space, is constant 0. We continue by creating the SPDE model through an R-INLA function inla.spde2.matern whose main arguments are used to pass the mesh object and to fix parameters α, τ and κ. Its arguments B.tau and B.kappa are matrices with one row for each mesh node. If only 1 row is given, it describes a model with stationary values of κ and τ, which will be duplicated internally for all nodes. For simulating a model with fixed parameters κ and τ, these matrices have only one column that must be specified as log κ or log τ respectively. Then, we extract the precision matrix Q Q Q of the resulting SPDE model and use it to create independent samples of ε(s,t), t = 1, . . . , 60 through the function inla.qsample. We fix the random seed for simulation through its seed=... argument. Finally, we manually create the first order AR model with coefficient 0.5.

:n.repl){x[,i]=a * x[,i-1]+sqrt(1-aˆ2) * x[,i]}
It remains to fix covariate values and to generate the time trend in the mean to add it to the centered Gauss-Markov space-time field. We fix an intercept -1 and the two covariates covar1 covar2. x y q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q FIGURE 2. For the observed data y y y to be used in estimation, we sample uniformly 50 sites among the triangulation nodes contained in [0, 1] 2 . By using R-INLA's methods inla.mesh.projector andinla.mesh.project to project a spatial field with known values for triangulation nodes onto a regular grid necessary for standard plotting methods, we further provide a plot of W (s, 19) and the observation sites. At t = 19, we observe the maximum value of the time trend, and we will later use R-INLA to do spatial prediction for t = 19. .mesh.projector(mesh.sim,dims=c(n.grid,n.grid)) image (grid$x,grid$y,inla.mesh.project(grid,field=x[,t.pred]),xlab="x", ylab="y",asp=1) points(sites,pch=19, cex=.5)

Soumis au

Figure 2 shows the trend component -1 + z 1 (t) + 0.5z 2 (t) (middle display) and the spatial field W (s, 19) at a fixed instant t = 19 with observation sites indicated (right display).

Fitting a complex space-time model with INLA

We now define and fit different candidate models to the above data y and its covariates covar1 and covar2. In Section 4.4, we will then explore tools for goodness-of-fit checks and model selection within the R-INLA package. We will first consider a model with prior structure similar to the simulated one, but we will also compare it to models where covariates are missing, where the Matérn shape parameter ν takes a different value or where the likelihood is not Poisson, but of the negative binomial type (with an additional hyperparameter for overdispersion).

First, let us define the triangulation mesh and the corresponding prior spatial SPDE model for estimation. For estimation, we must be careful about the dimension of the latent model to minimize memory requirements and high-dimensional matrix calculations. Therefore, will use a mesh with a lower resolution than for simulation, which may slightly increase the approximation error with respect to the stationary Matérn correlation structure. It is often reasonable to use observation sites as initial nodes of the triangulation and to refine it by adding nodes where necessary, or by removing nodes too close together which could be source for numerical instabilitities. R-INLA implements standard methods from the finite element literature and offers a conveniently parametrized interface to produce numerically stable and moderately dimensioned triangulations. R-INLA's prior parametrization of κ and τ is a bit technical; it essentially assumes that log(τ

) = b τ,0 + θ 1 b τ,1 + θ 2 b τ,2 + ... and log(κ) = b κ,0 + θ 1 b κ,1 + θ 2 b κ,2 + ...,
where the values b i correspond to the columns of B.tau and B.kappa. The first column b 0 is a fixed offset that must always be specified (even if it is 0), whereas the following columns correspond to hyperparameters θ i that are estimated when b i = 0 for τ or κ. For instance, specifying B.tau=B.kappa=matrix(c(0,1),nrow=1) would lead to a model where κ = τ = exp(θ 1 ). In our model, we fix the offset b 0 = 0 and we estimate two hyperparameters, one corresponding to log(τ), the other to log(κ). This can be seen as the standard prior specification of the SPDE model in R-INLA. We fix the "correct" simulated value of α in the SPDE model. mesh=inla.mesh.2d(sites,offset=c(-0.125,-.25),cutoff=0.075,min.angle=21, max.edge=c(.1,.25)) plot(mesh) points(sites,col="red",pch=19,cex=.5) spde=inla.spde2.matern(mesh=mesh,alpha=alpha,B.tau=matrix(c(0,1,0),nrow=1), B.kappa=matrix(c(0,0,1),nrow=1))

The mesh counts 294 nodes. Notice that further arguments of the inla.spde2.matern(...) function can be set to modify default priors, impose "integrate-to-zero" constraints, etc. We now run a first estimation by considering a model with Poisson likelihood and prior structure of the latent Gaussian field corresponding to the model that we simulated to generate the data. Naturally, this model should provide a good fit and we will later compare it to a number of alternative models. The observation matrix A A A links observations y y y to the latent variables x x x through A A Ax x x = y y y and must therefore be of dimension (50 × 60)×(number of latent variables) for our model. Since construction of this matrix and certain preprocessing steps before estimation like the removal of duplicate rows is rather complicated for complex models involving the spatial SPDE solution, R-INLA has helper methods that allow constructing this matrix and keeping track of latent variable indices more easily. In the following, the inla.spde.make.index command creates an index named "spatial", i.e., a data frame with a vector spatial (an index to match latent variables and triangulation nodes), a vector spatial.group (an index that indicates the membership of a latent variable in a "group" , here given as the instant t ∈ {1, . . . , 60}), and a vector spatial.repl (an index that indicates the group membership when groups are i.i.d. replicates). In our case, all values of spatial.repl are 1 since there is no structure of i.i.d. blocks in our space-time model due to the non-zero autoregression coefficient. idx.spatial=inla.spde.make.index("spatial",n.spde=spde$n.spde,n.group=n.repl) A.obs=inla.spde.make.A(mesh,loc=sites,index=rep(1:nrow(sites),n.repl), group=rep(1:n.repl,each=nrow(sites))) stack.obs=inla.stack(data=list(y=y),A=list(A.obs,1),effects=list(idx.spatial, data.frame(intercept=1,covar1=rep(covar1,each=n.sites), covar2=rep(covar2,each=n.sites))),tag="obs")

In practice, we may want to use a fitted model for prediction at unobserved sites. A natural way to achieve prediction in the Bayesian framework of INLA is to add the prediction points to the data by considering the associated observations as missing data. To illustrate R-INLA's facilities for this approach, we here consider prediction at instant t = 19 over a regular spatial grid covering [0, 1] 2 . Therefore, we first create a separate observation matrix and a stack for the prediction points with missing data, and we will then use R-INLA's join-mecanism that allows regrouping several groups of predictors η η η k through their observation matrices A A A k . The inla.stack.join(stack1,stack2,...)-function creates a single model corresponding to a structure (A A A 1 , A A A 2 , ...) x x x = (η η η 1 , η η η 2 , ...) , where information relative to A A A k x x x = η η η k is regrouped in stackk for k = 1, 2, . . . We now run the estimation with the inla(...) function. Its syntax ressembles that of R's glm(...)-function for generalized linear models, although with a variety of extensions and additional arguments. We need a model formula given in the usual R notation. For better handling of the intercept term, it is often preferable to make it appear explicity (form=y -1+intercept+...), such that it can later be directly included into the latent space(-time) model. Fixed effects (i.e., the intercept and covariates whose linear regression coefficients are estimated) are added to the formula in the usual way, whereas random effects are added with the f(...) function. In our model, the approximate SPDE solution is a random effect. The first argument of f is the name of the "covariate" associated to the random effect. Having created an index with name spatial beforehand, we now have a covariate spatial in the data set that indicates the triangulation node index of the spatial SPDE model (repeated 60 times since the spatial model is duplicated for each observation time). For prior SPDE models, we further specify the argument model=spde in f, where spde is the R object already created for the SPDE prior model. The SPDE model is of purely spatial nature whereas we have observations in space and time, such that we can use the group-functionality of R-INLA to define the type of dependence between the 60 groups of spatially indexed Gaussian variables. Corresponding to the simulated model, we here use an AR(1)-group model that models site-wise first-order autoregression of variables over time. Since we have created the index spatial, we can specify the argument group=spatial.group to indicate group membership of the covariates, and we have to set control.group=list(model="ar1") for the AR(1)-model between groups.

Data must be passed to the fitting function inla(...) as a data.frame or list, and the inla.stack.data-function allows convenient extraction of the preprocessed data object from the stack. Further control arguments to inla(...) can be specified explicitly through R's usual control. =list(...) syntax, which allows overriding the default control arguments. Here, should be replaced by one of inla (for controling INLA-related details like the reordering scheme used for making the precision matrix as diagonal as possible), compute (for specifying which quantities should be calculated, e.g. goodness-of-fit and model selection criteria), predictor (for specifying the observation matrix A A A if there is one, and for indi- cating which posterior quantities should be calculated for the predictor vector η η η), family (for modifying the default prior of likelihood hyperparameters), amongst others. In the following, we fix the METIS-reordering strategy in control.inla to avoid instabilities that have appeared during a preliminary fitting approach with the default reordering scheme. In control.predictor, we pass the observation matrix A A A that can be extracted from the stack via inla.stack.A(stack), and we advise the program to calculate (discretized) posterior densities for the predictor variables η i (compute=T). In control.compute, we here demand the calculation of CPO-values (cross-validated predictive measures, see [START_REF] Held | Posterior and cross-validatory predictive checks: a comparison of MCMC and INLA[END_REF] for a comparison of Markov chain Monte Carlo and INLA), the marginal likelihood π(y y y), the Deviance Information Criterion (DIC) and the Watanabe-Akaike Information Criterion (WAIC, [START_REF] Watanabe | Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory[END_REF]; [START_REF] Gelman | Understanding predictive information criteria for bayesian models[END_REF]), where the latter can be considered as a Bayesian variant of the common AIC. The CPO-related values are the density and cdf of the posterior distribution π(y i | y y y -i ), evaluated at the observed y i . With INLA, these cross-validation quantities can be calculated quickly without explicitly reestimating the model, see [START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF] for details. However, certain theoretical assumptions might be violated such that some or all of these CPO-related values are not trustworthy, which is then indicated in the inla-output in fit$cpo$failure. In such a cases, the inla.cpo(...)-function can be used for "manual" reestimation of the cross-validated model for the concerned data points y i . We now construct the inla(...)-call. For later use, we here store the data and parameters of this first model in an object mod1: data=inla.stack.data(stack,spde=spde) form=y˜-1+intercept+covar1+covar2+f(spatial,model=spde, group=spatial.group,control.group=list(model="ar1")) cc=list(cpo=T,dic=T,mlik=T,waic=T) cp=list(A=inla.stack.A(stack),compute=T) ci=list(reordering="metis") mod1=list (stack=stack,data=data,A.pred=A.pred,A.obs=A.obs, idx.spatial=idx.spatial,spde=spde,form=form,cp=cp,ci=ci,cc=cc) fit=inla(form,family="poisson",data=data,control.compute=cc, control.predictor=cp,control.inla=ci) Since our likelihood is not Gaussian (but Poisson) and since the latent model is quite complex, the inla run takes some time (around 50 minutes on a state-of-the art 4 core machine), and memory requirements are rather high. We remark that inla(...) has a num.threads argument which allows fixing the maximum number of computation threads used by R-INLA. By default, R-INLA takes control over all available cores of the machine for parallel execution, which can lead to problems in terms of too high memory requirements since each thread occupies a certain amount of memory. As is usual in R, we can now call summary(fit) to obtain principal results of the fit. Part of its output is as follows: Cdf values are not very far from being uniform, which indicates that there are no strong systematic biases in posterior predictions made with the model. A certain proportion of the failure indicator values are above zero and some are even 1, meaning that there are some cdf values that must be interpreted with caution.

Since we have done prediction for t = 19, we now extract and visualize the marginal posterior mean and standard deviation over the prediction grid. To get the index of predicted variables, we can apply the inla.stack.index function to the stack object by indicating the prediction sub-stack through the argument tag="pred". The following code visualizes the predictions ηi and the originally simulated Gaussian values on the prediction grid. The x=xgrid,y=xgrid,eta.mean,asp=1,xlab="x",ylab="y",main="posterior mean") proj=inla.mesh.projector(mesh.sim,xlim=c(0,1),ylim=c(0,1),dims=c(n.grid,n.grid)) image(grid$x,grid$y,inla.mesh.project(proj,field=x[,t.pred]),xlab="x", ylab="y",xlim=c(0,1),ylim=c(0,1), asp=1,main="original") points(sites,pch=19,cex=1) eta2.mean=unlist(lapply(eta.marginals,function(x) inla.emarginal("ˆ",x,2))) eta2.mean=matrix(eta2.mean,n.grid,n.grid) eta.sd=sqrt(eta2.mean-eta.meanˆ2) image (x=xgrid,y=xgrid,eta.sd,asp=1,xlab="x",ylab="y",main="posterior sd") points(sites,pch=19,cex=1)

Figure 4 shows the resulting displays. As expected, uncertainty is lower close to observation sites. The prediction captures the spatial variation of the actual data, and a deeper analysis shows that the predicted surface is smoother than the original values: since the Gaussian prior on the predictor is centered at 0, the spatial variation in posterior predictions is naturally less strong owing to the Bayesian approach.

Other candidate models

In Section 4.2, we have used our knowledge about the simulated data model to construct the Bayesian model that should be the most appropriate. For comparison and to illustrate R-INLA's functionality for other types of models, we here propose to fit some alternative candidate models.

For Model 2, we drop the temporal autoregression and consider the spatial blocks of SPDE variables as independent in the prior. This necessitates some (slight) adaptations in the code since Finally, Model 6 uses not the Poisson likelihood but the negative binomial one that has an additional overdispersion parameter θ disp , where the variance of the negative binomial distribution is µ + µ 2 /θ disp and µ is its mean. Notice that the Poisson distribution arises in the limit for θ disp = ∞. The internal parametrization of θ disp considers log θ disp as a hyperparameter. We here fix a relatively high initial value log(10) and use a log-gamma prior with shape 10 and rate 1, such that the prior expectation of θ is 10 and yields a prior likelihood model relatively close to the Poisson one: cf=list(hyper=list(list(theta=list(initial=log(10),prior="loggamma", param=c(10,1))))) fit=inla(mod1$form, family="nbinomial",data=mod1$data,control.compute=mod1$cc, control.family=cf,control.predictor=mod1$cp, control.inla=mod1$ci)

Analyzing fitted models

We now compare the 6 fitted models. For the purely temporal model 3, the following code plots a posterior mean estimate of the fitted temporal curve and the simulated data, see Figure 5:

plot(1:n.repl,fit$summary.fixed[1,1]+fit$summary.random$covar1$mean,type="l", xlab="t",ylab="time trend",lwd=2) lines(1:n.repl,xtrend,col="blue",lwd=2)

We see that neglecting the spatial variation in data and using a relatively simple, purely temporal model here permits to reconstruct quite accurately the simulated temporal trend. Finally, we can compare the information criteria DIC and WAIC, marginal likelihoods, estimates of spatial range, variance and temporal autocorrelation over the candidate models, see Table 1. The marginal likelihood is relatively close for all space-time models with an explicit spatial structure and temporal autoregression (Models 1,4,5,6), and has its by far lowest value for the purely temporal model 3. Model 1 whose prior structure is closest to the simulated model turns out best in terms of DIC, but has slightly higher WAIC values than models 4 and 5 characterized by a different fixed shape parameter ν in the SPDE solution, leading to less smooth sample paths in the spatial prior random field. Notice however that differences in the estimated WAIC values between models 1,4 and 5 are relatively small such that they should be interpreted with caution. Model 6 with the negative binomial likelihood but with the same latent Gaussian prior model as Model 1 has rather high DIC and WAIC values, maybe due to approximations and computations that are less accurate for this fitted model. Concerning range and variance parameters of the fitted spatial models, we find that they are close to simulated values in all cases except Model 5, where the different prior shape parameter ν = 0 in the SPDE seems having perturbed the calculations of the dependence structure and the variance. Throughout, the posterior distribution of the autoregression coefficient a (if estimated) is concentrated around the true value 0.5.

Discussion

We have illustrated the theory and practice of Integrated Nested Laplace Approximation, implemented in the powerful R-INLA software library, which enables fast and accurate inference of complex Bayesian models. The dynamic behavior and dependence structure in the models covered by R-INLA is primarily governed by a latent Gaussian process for the mean of the univariate data distribution. In view of the near endless range of models that are available, one should perhaps also sound a note of caution since users may be mislead to construct overly complicated models in practice.

R-INLA has mechanisms to manage several different likelihoods in the same model, to use the same latent variables in different parts of the latent model ("copying") or to deal with noninformative missing data. Although the INLA approach does not directly provide the posterior dependence structure between predictors η i or between the data distributions for different y i , one could assume a Gaussian copula model with the Gaussian dependence given by the precision matrix Q Q Q * of the Gaussian approximation (11) to obtain a practically useful approximation of the posterior dependence. Moreover, R-INLA has inbuilt functionality for calculating posterior densities of user-defined linear combinations of the latent variables, i.e., for latent Gaussian variables other than the predictors η i or the components of the latent model; see the lincombargument of inla(...). New users should look around at www.r-inla.org, the main hub for staying informed about new INLA-related developments, for finding implementation details of R-INLA and for getting advice on specific problems via its discussion forum.

The SPDE approach, providing flexible spatial Gauss-Markov models, is of interest in its own far beyond the INLA framework where Markovian structures lead to fast high-dimensional matrix computations. Multivariate extensions [START_REF] Hu | Multivariate Gaussian random fields using systems of stochastic partial differential equations[END_REF] and certain nested version of SPDEs [START_REF] Bolin | Spatial models generated by nested stochastic partial differential equations, with an application to global ozone mapping[END_REF] have already been proposed in the literature, although they are not yet available in R-INLA. Further modeling extensions in terms of data likelihoods and latent models can be expected in the near future. In particular, the construction of more complex and realistic Gauss-Markov space-time dependence structures based on the SPDE approach, capable to model effects like the nonseparability of space and time would be another big step forward.

FIGURE 1 .

 1 FIGURE 1. Illustration taken from[START_REF] Blangiardo | Spatial and Spatio-temporal Bayesian Models with R-INLA[END_REF]. Example of a Gaussian surface (left) and its finite element approximation (right). On the right display, we can further see the triangulation and one of the finite element pyramid-shaped basis functions.

  sqrt(pi) * kappa * sigma0)

Constrained

  Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016

  and the "empirical range" where a correlation of approximately

Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016

  Simulated latent Gaussian model. Left: triangulation used for the approximate SPDE simulation. Middle: time trend. Right: simulation for t = 19, with observation sites marked by black dots.

	covar1=1:n.repl/n.repl
	covar2=as.numeric(arima.sim(n=n.repl,model=list(ma=rep(1,5))))
	xtrend=-1+covar1+0.5 * covar2
	x=t(t(x)+xtrend)
	plot(xtrend,type="l",xlab="time",ylab="trend",lwd=2)

  Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016

  Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016

  Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016 Left: histogram of Internally cross-validated cdf values for π(y i | y y y -i ). Right: histogram of the failure indicator of trustworthiness of calculated cdf values.inla.mesh.projector function allows switching between the finite element representation and a regular grid by calculating the Gauss-Markov finite element approximation value to the SPDE for the grid positions by using the "pyramid" basis functions to interpolate between triangulation nodes. We use the inla.emarginal-function for calculating posterior expectations, where the standard deviation is calculated as EX 2 -(EX) 2 in the following code:

	number	0 50 150 250					number	0 500 1000 1500
		0.0	0.2	0.4	0.6	0.8	1.0	0.0	0.2	0.4	0.6	0.8	1.0
				PIT value					failure indicator
	FIGURE 3.

idx.pred=inla.stack.index(stack, tag="pred")$data eta.marginals=fit$marginals.linear.predictor[idx.pred] eta.mean=unlist(lapply(eta.marginals,function(x) inla.emarginal(identity,x))) eta.mean=matrix(eta.mean,n.grid,n.grid) image(

TABLE 1 .

 1 Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016 Purely temporal random walk model 3. Posterior mean curve (black) and simulated curve (blue). Comparison of fitted models: DIC, WAIC, marginal likelihood mlik, spatial range, variance of spatial model, autoregression coefficient a. All values are rounded to three significant digits and to at most two decimals. Estimates are posterior means and 95% credible intervals are in parentheses. Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016

				3
				2
				1
			time trend	0
				-1
				-2
				-3
				0	10	20	30	40	50	60
							t
	FIGURE 5. Model DIC WAIC		mlik			range	variance	a
	1	8540	8520 -4750 0.29(0.25;0.33) 0.95(0.84;1.08) 0.49(0.43;0.55)
	2	8660	8630 -4830 0.30(0.26;0.34) 0.92(0.82;1.03)	-
	3	9530	9270 -5860			-	-	-
	4	8550	8510 -4750 0.32(0.26;0.38) 1.03(0.91;1.16) 0.49(0.42;0.55)
	5	8550	8500 -4760 0.59(0.43;0.79) 0.11(0.09;0.14) 0.48(0.42;0.55)
	6	8890	8810 -4750 0.34(0.29;0.40) 0.82(0.70;0.96) 0.53(0.46;0.60)

Soumis au Journal de la Société Franc ¸aise de Statistique File: INLASFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 26, 2016

T. Opitz x y q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 posterior sd x y q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q FIGURE 4. Predictions of the latent Gaussian field for t = 19. Left: original simulated data. Middle: posterior mean. Right: posterior standard deviation.

there is no more group model with dependence between blocks, but we now have replicates, i.e., independent blocks. Model 3 focuses on the temporal trend, neglects spatial variation and supposes that covariates are not known. Here we use a temporal first-order random walk as prior model. We set a rather high initial prior value for the precision of the random walk innovations (corresponding to a standard deviation of 0.01). A sum-to-zero constraint is added (constr=T) for better identifiability.

data3=data.frame(y=y,covar1=rep(covar1,each=n.sites)) form=y˜f(covar1,model="rw1",hyper=list(prec=list(initial=log(1/.01ˆ2), fixed=F)),constr=T) cp=list(compute=T) fit=inla (form,family="poisson",data=data3,control.predictor=cp) Another interesting model could be to obtained from combining Models 2 and 3, i.e., using a random walk in time and a spatial SPDE model. Such a model is relatively complex if not overly complicated and its estimation is computationally demanding, so we do not consider it here. In Models 4 and 5, we specify a shape parameter of the Gauss-Markov Matérn model in the prior that is different from the simulated model, using either ν = 0.5 (exponential model) or ν = 0 (not a proper Matérn model, but still a valid covariance model).


Martin Petit

Eric Prada

Valérie Sauvant

Valérie Sauvant-Moynot

Development of an empirical aging model for Li-ion batteries and application to assess the impact of V2G strategies on battery li/ etime

Development of an empirical aging model for Li-ion batteries and application to assess the impact of V2G

strategies on battery lifetime

Introduction

ln order to reduce fossil energy dependence and environmental impact of vehicles, new regulations encourage car manufacturers towards vehicle electrification. This electrification raises the issue of energy storage. Li-ion batteries are one of the most promising solutions to store the energy needed for highly electrified vehicles: hybrid electrical vehicles (HEV), plug-in hybrid vehicles (PHEV) or full electric vehicles (EV). However, energetic and power performance of Li-ion batteries is known to decrease during their service lifetime (1).

Although the cost of battery packs is decreasing [2], they comprise a significant proportion of market vehicle price, which means car manufacturers have to provide a guaranteed longevity ensuring a safe and sufficient battery performance throughout the vehicle service lifetime. Therefore, battery aging studies are of particular interest in order to optimize battery systems in terms of size and cooling requirements, and also to devise management strategies well suited along the vehicle lifetime. Providing a minimum level of performance is also mandatory to ensure financial viability for electrified vehicles. For such concerns, modeling is essential to provide car designers with information and reduce development costs.

White the main microscopie aging phenomena of Li-ion cells have been identified [1,3) and described, for instance Li plating at low temperature (4), loss of active material at elevated temperature (5) or solid electrolyte interphase (SEI) growth at the negative electrode/electrolyte interface [6,7), aging of Li-ion batteries remains difficult to predict [8,9]. ln most aging models, considering that SEI growth is the main cause of aging for Li-ion batteries with graphite negative electrode leads to a quadratic evolution shape of the capacity loss in the beginning of the battery lite as put forward by Spotnitz (10).

To find the main impact factors, research efforts have been dedicated to model macroscopic Liion cell capacity loss and impedance increase [START_REF] Wang | 1-1[END_REF]12). The main aging factors are the state of charge, SoC, the depth of discharge during cycle, DoD, the temperature, T, and current, /. Another key impact factor is the battery usage (13) which can be described in two contributions to aging : Calendar aging (/ = O): the battery is stored without being used so there is no current through the battery Cycle aging (/ -:/: 0): the battery is either charged or discharged. lt has been shown that the two types of aging will lead to different battery behaviors with usually a higher capacity loss and higher power loss during cycle aging (13].

ln order to model aging phenomena, two different approaches have been put forward in the literature, which are physical modeling [START_REF] Newman | Electrochemical Systems[END_REF][START_REF] Prada | [END_REF](16) and empirical modeling (10,17]. Usually empirical models will address either calendar aging (10,13] or cycle aging [START_REF] Wang | 1-1[END_REF]. Sorne other modeling works include both cycle and calendar aging behaviors [18][19][20] but the proposed model despite describing cycle aging does not take into account the charging rate on aging. The resulting model assesses the capacity loss as a function of time and constant stress factors which has little application in use. lt is not convenient to use this kind of modeling for system simulation where stress factors are dynamically varying for instance due to environmental conditions {temperature variations due to seasons and day and night successions) or due to vehicle use {storage Soc depending on the duty cycle). ln this paper, an empirical capacity fade model suitable for system simulation application is described. This model has been calibrated for two li-ion technologies and validated against experimental data and is used to assess the usefulness of different charging strategies in order to preserve the batteries' performance throughout their lifetime.

Mode/ development

The capacity fade model is based on 2 contributions leading to different capacity lasses in Ah:

Qf 0 a} 5 due to calendar aging and Qf:S~ due to cycle aging.

2.1.

Calendar aging

As often discussed the two main stress factors for calendar aging are temperature and State of Charge. Usually [START_REF] Wang | 1-1[END_REF]21) an empirical aging law is given as below (1) ln this expression Bcal is a pre-exponential factor depending on SOC, expressed in :h 1 , Eacal is s ca the activation energy, expressed in J.mor 1 , which evaluates the dependency of calendar aging on temperature T, expressed in K, and Zcal is a dimensionless constant. Considering a capacity loss phenomenon linked to SEI growth and diffusion limitations, this exponent should be a round O.S. This expression is well suited to rapidly evaluate the capacity lost during long term storage but it is not suited when many different operating conditions are happening. As a consequence, this expression has been differentiated against time in order to evaluate infinitesimal variation of capacity due to calendar aging:

1-1- dQcal E ( Qcal ) Zcal ~ _ B (SOC) (-acal) loss -Zcal cal exp E dt RT Bca1(SOC) exp (-~Tal) (2)
This expression can then be linked to other capacity loss phenomena like cycle aging.

Cycle aging

For cycle aging a similar approach based on Wang et al. [START_REF] Wang | 1-1[END_REF]

dQ~:S~ III (-Eacyc + alll) ( Q~:S~ )'-,~, (4)
~ = 3600 ZcycBcyc(/) exp RT (-Eacyc + al/I) Bcyc (/) exp
RT ln this expression, the Ah throughput term, once differentiated, is proportional to the current in the cell.

2.3.

Switch between calendar and cycle aging

Once the influence on aging has been expressed for both behaviors, it is necessary to find a way to decide what kind of aging is occurring based on the actual operating conditions. The assumption here is that when the modeling of cycle aging has been performed, calendar aging occurring during cycle mode is already taken into account. As a consequence, we need to find a condition describing the switch between calendar and cycle mode.

First, taking as a major phenomenon of capacity loss the SEI formation at the negative electrode, it is known that the parasitic reaction of SEI formation is accelerated when the negative electrode potential is lower (i.e. during charge or when the SOC is high). As a consequence it has been assumed that cycle aging only occurs when the battery is in charge and the current is above a given threshold Icyc since higher charge current will lead to lower potential at the negative electrode. Furthermore, in order to assess the fact that diffusion in the cell tends to increase the overpotential when the charge is longer and then increases the rate of the SEI formation reaction, the current in d~f~ expression ( 4) is replaced by a filtered current Ï. As a consequence the total capacity loss is computed as follows:

{ dQcal ~ ~Ï<l dQ1oss = dt cyc dt dQcyc ~~Ï>l dt cyc ( 5
)
This given threshold will depend on the battery ability to tackle high charging rates. For example high power batteries should have a high threshold around lC current whereas high energy batteries designed for lower charging rates should have a lower threshold.

2.4.

Coupling with an electrothermal mode

Experimental validation

Calibration data

ln order to calibrate the empirical model presented in this paper, extensive experimental data are needed in order to get a wide range of operating conditions representative of the battery use.

Experimental campaigns have been carried out by research groups in order to assess bath cycle and calendar behavior of Li-ion batteries [12,18,[20][21][22][23]. Two sets of calibration data are indeed necessary in order to take into account bath types of aging. ln these sets, experiments highlighting every stress factor are needed. As a consequence a significant amount of experimental data was needed and was collected in the literature for calendar aging [24) as well as for cycle aging [START_REF] Wang | 1-1[END_REF][START_REF] Gyan | Oil and Gas Science and Technology[END_REF].

Thanks to these data, it has been possible to calibrate the capacity loss model on two different battery technologies:

A123s 2.3 Ah LFP/C cell Saft VLGP 7 Ah NCA/C cell
The parameters used for these technologies are indicated in Table 1. ln the following section we will discuss in detail the validation for the A123s model.

Validation of empirical model

Normalized profiles

Extensive experimental studies from Savoye [START_REF] Savoye | Impact des impulsions périodiques de courant sur la performance et la durée de vie des accumulateu rs lithium-ion et conséquences de leur mise en oeuvre dans une application transport[END_REF] have subjected an A123s 2.3 Ah cell to various cycle conditions in order to evaluate the cycling impact on Li-ion batteries aging. The aim was to discuss the aging impact in an automotive application and especially the impact of short current pulses on aging. The current profiles represented in Figure 2 are looped at 25°C and 45°C over more than 4000 cycles. These current profiles interestingly caver various cycle profiles with highly dynamic pulses at several different rates (ISO RMS and ISO MEAN), moderately dynamic pulses (DILATED) and slow solicitations (REFERENCE or 0.8 REFERENCE), as well as mixed sa licitations (AT REST).

The comparison between experimental and simulation results using the empirical aging model described above is presented in Figure 3. lt shows a good agreement between experimental results

and simulations with less than 5% error for capacity loss lower than 20% except for the ISO MEAN condition where the model overestimates the capacity loss by 14%. These results seem encouraging for using this model in automotive applications.

3.2.2.

Road profile

ln order to further validate the empirical model on actual data representative of vehicle use, aging measurements have been performed at cell level on IFPEN test benches [START_REF] Bernard | Batteries: The International Energy & Power Supply conference and Exhibition[END_REF]. Two different load profiles representative respectively of EV and HEV were looped. For the EV load profile (Figure 4-a), the program consisted of 4 charge depleting cycles followed by a constant current charge in order to get back to the initial SOC. The HEV load profile (Figure 4-b) consists of a charge-sustaining current profile which maintains the SOC around 60%. The maximum current involved in the HEV profile is a round 50 A whereas it remains lower than 20 A in the EV profile.

ln Figure 5, it is possible to verify that the model is able to predict accurately capacity loss for bath aging patterns (EV in Figure 5-a and HEV in Figure 5-b). lt is then possible to use this model in automotive application.

Applications for lifetime prediction

This capacity loss empirical model has been implemented in LMS lmagine.Lab Amesim™ dedicated to system simulation. lt can then be used in realistic operating conditions in order to adapt the charging strategies taking into consideration the long term issue of aging. ln this scope we focused on a study that will discuss the impact of the recharge strategy of an electric vehicle fitted with a LFP-C battery or a NCA-C battery.

Case studies

Day to day use of an electric vehicle can be simplified with a common pattern comprising charging phases and trips from home to work and back. The base case which is considered in this study is the following:

4 home-work trips lasting 20 minutes at SAM, 12PM, lPM and GPM.

ln addition to these trips several charging and discharging phases are added. The discharging phase represents the vehicle connected to grid sections where the energy stored in the vehicle is used to relieve the grid during high demand periods. The constant charges and discharges are performed at a 40 A rate and stopped whenever the SOC reaches 10% or 90%. As a consequence 5 test cases scenarios have been defined (Figure 6):

The reference case consists in a classical daily use with a simple recharge when the driver cornes back home.

The "just in time" case consists in charging the car just before the first trip in the morning.

During the "charge when you can" case the driver tries to charge the car whenever he has time to do so (that is after every trip except when he gets back home for lunch).

The "strong V2G" scenario consists in a discharge of the car when there is time (at 9AM, 2PM

and 19PM) and a recharge when the battery is empty before leaving (GAM, llAM, 4PM). As a consequence, the battery is mostly stored at a low Soc (10%).

The "light V2G" scenario is only considering a discharge of the battery when the driver cornes back home at the end of the day.

These test case scenarios deal with two main stress factors regarding aging, which are current inside the cells and the SOC. Another stress factor of aging to be accounted for is the temperature.

Figure 7 illustrates the temperature variations seen by the battery during a year in this case study.

There is in tact a double temperature variation. External temperature is submitted to a daily variation due to days and nights but also to an annual variation due to seasons. Moreover, the use of the battery leads to a temperature increase due to heat losses. As can be seen in the prediction of battery daily temperature variations some small peaks mark a battery self-heating of 1 or 2°C under the hypothesis of an efficient thermal exchange leading to an exchange coefficient of 50 W/m 2 /K.

Results

Simulations have been performed on both LFP/C and NCA/C packs whose properties are given in Table 2. The aim was to build a battery pack with the sa me capacity. However due to the differences between the nominal voltage and inherent properties of the cells such as internai resistance, the pack power and energy are different.

The simulation evaluates the aging of the batteries du ring one year use considering that the profile patterns described in Figure 6 occur every day. ln Figure 8 are represented the capacity losses for both technologies for every test case.

Due to the sizing differences between both packs, we will not focus on the comparison of the raw performances of both technologies but rather on the differences that can be seen between the test case behaviors.

For both cells, the best scenario is the "just intime" case. ln this scenario, the cell remains at low SOC during the night, reducing the calendar aging effect compared to the reference case where it is stored at 90% SOC. lt is even more visible in the "charge when you can" scenario for the LFP/C based battery where the cell is always stored at 100% SOC increasing capacity loss during storage. Yet the "strong V2G" scenario where the battery is always stored at 10% SOC is not so good compared w ith the reference case for both technologies. ln this case the advantage from storage at low SOC is outweighed due to the need to recharge the battery to 100% SOC before every use. This recharge causes a higher cycle aging capacity loss and th us a premature aging.

However, the behavior seen in the case of a LFP/C based battery is different from the one encountered with a NCA/C based battery where the aging of the "charge when you can" scenario is slightly bigger than the reference case but stays in the same order of magnitude. Moreover it clearly appears that in the case of the NCA/C based pack, the V2G scenarios lead to a premature aging of the battery (4% loss for the strong V2G and 2.8% for light V2G) and are even worse than the "charge when you can" scenario (2.7% loss) whereas this last scenario is the worst one in the LFP/C based battery case (7.6% loss) and the light V2G is the second best scenario for the LFP/C battery on par with the "just in time" case (a round 3.5% loss). lt seems that the NCA/C based battery studied here is more sensitive to cycle aging than calendar aging. Then it will perform worse when submitted to frequent charges (eg. during V2G and "charge when you can" scenarios).

As a consequence, both technologies exhibit clearly different behaviors when submitted to the sa me test cases. Based on these simulations, advice for vehicle use and especially recommendation on V2G applications will be different depending on the technology used to build the battery pack.

Conclusion

An empirical capacity aging model for Li-ion batteries has been developed. lt is coupled with an electrothermal model and accounts for both types of aging encountered during Li-ion battery life: calendar aging and cycle aging. Thanks to an extensive experimental campaign and data mining it has been possible to calibrate and validate the model for two different Li-ion technologies. Once validated, the model has been integrated in a system simulation software application in order to assess the behavior of the cells when submitted to various test case scenarios. These simulations showed some major differences between both technologies, the NCA/C based battery proved to be more sensitive to cycle aging compared to the LFP/C cell leading to a premature aging when submitted to frequent charges of the V2G scenarios.

This approach proved to be appropriate in order to evaluate the charging strategies that could be applied to an electrical vehicle in order to reduce its aging and in doing so increase its lifetime and cost efficiency. lt is then able to accurately model the specific behaviors of both cells used in this study leading to recommend a light V2G strategy for the LFP/C based cell when it should rather be avoided in the case of the NCA/C based battery.

This study only focuses on the capacity loss due to aging. However, it is known that batteries also exhibit resistance increase. This increase has been neglected here since in the case of both cells used this resistance increase remained low even at the end of life. However, for other cells this may not be the case. As a consequence further work should be done in order to apply this approach for power loss.

  / of the battery behavior ln order to get a global model of the electrothermal behavior of a Li-ion cell, it is necessary to couple this aging model with an electrothermal model. ln Figure 1 are represented the inputs and outputs of both models as well as the coupling involved. The electrothermal model evaluates the state of charge (SOC), the voltage (Veen) and the temperature (Tceu) of the battery based on the input current (1) and external temperature (Ta). lt provides the empirical aging model with the SOC and the Tcen allowing it to evaluate the residual capacity Cres • This residual capacity is then used in the electrothermal model to evaluate the SOC as follows:100 f t I SOC = SOCo + 3600 o Cres(t) dt(6)

1 Figure 1 :Figure 2 :

 112 Figure 1: Representation of the coupling between an electrothermal model and the empirical aging model

Figure 4 :

 4 Figure 4: Current profile used for EV and HEV validation

Figure 5 :

 5 Figure 5: Comparison between experimental data and simulation results

Figure 6 :

 6 Figure 6: Current and SOC profiles during test case scenarios

NCNCFigure 8 :

 8 Figure 8: Capacity loss for both battery technologies after 1 year cycle according to every scenario

  Bcyc is a pre exponential factor in Ah 1 -zcyc which depends on current, Eacyc Ah stands for Ah throughput, that is the amount of charge sent into the cell. As it has been done for calendar aging, this expression has been differentiated against time:

	cyc Qloss = Bcyc exp	(-Eacyc + al/I) RT	AhZcyc	(3)
	ln this expression, is an activation energy for cycle aging expressed in J.mor 1 , ais a coefficient for aging acceleration
	due to current expressed in J.mor 1 A1, Zcyc is an exponent constant that should be around 0.5 for
	diffusion limited process. Finally,
			is used. ln this study the capacity
	loss due to aging is influenced by two major stress factors: current and temperature . The formulation
	encountered is then:

Table 1 :

 1 Parameters of the empirical aging model

		A123s 2,3Ah			SaftVL6P
	bO c •a.o "' Eaca1(SOC) Bca1(SOC) ... "' "tJ c Qj Zca1(SOC) "'iii u	30% 734063 30% 73369J/mol 30% 0.943	65% 675334 65% 69804J/mol 65% 0.900	100% 218211 100% 56937J/mol 100% 0.683	30% 2.78 • 10 5 1 3.80 • 10 5 1 5.43 • 10 5 1 65% l 100% 52862 0.52
	Bcyc(/)	lA	l4A	1 12 A	1 20A	130
	bO c ïio "' bO Zcvc Eacvc .5 ü acvc ('; fcvc	3.16. 10 3 1 2.17. 10 4 1 1.29 • 10 4 1 1.55 • 10 4 31700J/mol 0.55 370.3 lC	18461J/ mol 0.4 32 lC
	Tfilt	605				Os

Table 2 :

 2 Battery pack properties for simulations

	Chemistry of battery cells	LFP/C	NCA/C
	Capacity of the battery	2.3Ah	2.3Ah
	Number of cells in series	100	100
	Number of parallel branches	10	10
	Power of the pack	384kW	189kW
	Energy of the Pack	7600Wh	8280Wh
	Nominal Voltage of the pack	320V	360V