Analysis, Simulation and Prediction of Multivariate Random Fields with Package RandomFields

Modeling of and inference on multivariate data that have been measured in space, such as temperature and pressure, are challenging tasks in environmental sciences, physics and materials science. We give an overview over and some background on modeling with cross-covariance models. The R package RandomFields supports the simulation, the parameter estimation and the prediction in particular for the linear model of coregionalization, the multivariate Mat´ern models, the delay model, and a spectrum of physically motivated vector valued models. An example on weather data is considered, illustrating the use of RandomFields for parameter estimation and prediction.


Introduction
Spatial data very frequently have more than one component.In meteorology, for instance, temperature, rainfall and pressure are measured at the same locations at predefined instances in time.Similarly, in order to estimate the spatial distribution of rainfall, data from rain gauges are combined with radar data.In economy, to name a further example, price developments of goods are spatio-temporal data consisting of various components at each instant of time and at each market.
Although these kinds of data are ubiquitous, multivariate spatial models are largely underdeveloped and multivariate space-time models are unknown, except for simple constructions.This paper deals with multivariate spatial data, only.We denote a random process on a subset of R d , d ≥ 1, with scalar values a univariate random field, whilst other authors denote such a process already multivariate if d ≥ 2. Here, a multivariate random field is an R m -valued random process on a subset of R d with m ≥ 2. Such fields are called multivariable in the package gstat (Pebesma 2004).
The analysis and the modeling of data stemming from multivariate random fields have been reported in the geostatistical literature since long, cf.Wackernagel (2003) and Chilès and Delfiner (1999), for example.Recently, the subject of multivariate geostatistics has been taken up and various new models have been suggested (Apanasovich and Genton 2010;Gneiting, Kleiber, and Schlather 2010;Bevilacqua, Daley, Porcu, and Schlather 2013).
While spBayes mainly focuses on spatial regression models, multivator allows for inference and prediction of multivariate data based on emulator techniques (Oakley 1999).The package gstat provides several tools for prediction and model fitting and also allows for multi-Gaussian sequential simulation.However, to the best of our knowledge, the package RandomFields (Schlather et al. 2014b) is by now unique in its variety of tools and models for multivariate spatial processes.
The paper is organized as follows.In Section 2, a summary of the theoretical background is given.Section 3 discusses principles for analyzing and simulating multivariate random fields.Section 4 offers an overview over multivariate and vector-valued covariance functions that are implemented in RandomFields.Section 5 reconsiders the data analysis by Gneiting et al. (2010).Section 6 concludes with a brief summary of other features of RandomFields that are already available or are intended to be implemented.

Background: cross-covariance functions and cross-variograms
Throughout the paper, a second order m-variate random field on T ⊂ R d denotes a collection of real-valued random vectors Z(x) = (Z 1 (x), . . ., Z m (x)) indexed by x ∈ T with existing second moments.For theoretical considerations it suffices to assume EZ j (x) = 0 for all x ∈ R d and j = 1, . . ., m.The cross-covariance function C : T × T → R m×m is defined by C jk (x, y) = COV(Z j (x), Z k (y)), x, y ∈ T, j, k = 1, . . ., m.
The function C is positive definite, i.e., C(x, y) = C (y, x) and for all n ∈ N, x, y, x 1 , . . ., x n ∈ T and a 1 , . . ., a n ∈ R m .The cross-variogram γ : T × T → R m×m , is defined by Clearly, stationarity implies intrinsic stationarity.Further relations between covariance functions and variograms exist, see Chilès and Delfiner (1999) and Schlather (2012), for instance, and the references therein.
In the univariate case, (weak) isotropy of Z is uniquely defined and means that the function C or γ is rotation invariant.In the multivariate case the notion of isotropy has different interpretations.Assume that the bivariate field Z models temperature and pressure.Then Z should be called isotropic if Assume now that Z models a vector field of directed quantities in R d , e.g., wind speed in R 2 .Then Z has necessarily m = d components and it should be called isotropic if Z(•) A similar definition holds for the variogram γ.We call the latter models isotropic vector-valued models whilst the former ones are referred to as isotropic multivariate models.

Simulation and inference
Several methods for the simulation of univariate random fields in R d exist, see Lantuéjoul (2002) for an overview.The simulation of univariate random fields based on Cholesky decomposition, see also SpatialTools, can be directly transferred to the multivariate case.As the Cholesky decomposition currently reaches its limits with a matrix of about size 10 4 × 10 4 , only about 2500 locations can be modeled in a bivariate setup and about 1000 locations in a trivariate setup.The circulant embedding method by Dietrich and Newsam (1993) and Wood and Chan (1994) is extremely powerful for simulating on a regular grid in R d when d is small and the random field is stationary.Their algorithms are implemented in various packages, e.g., dvfBm (Coeurjolly 2009), fields, and stpp (Gabriel, Rowlingson, and Diggle 2013).
The method has been extended by Dietrich and Newsam (1996) and Chan and Wood (1999) to the multivariate case.In RandomFields, the Cholesky decomposition and the circulant embedding method suggested by Dietrich and Newsam (1996) are implemented also in the multivariate case.Furthermore, an extended version of Matheron's turning bands method (Matheron 1973) for the multivariate case is implemented.
Various methods exist to estimate the parameters of an intrinsically stationary random field.
Least squares methods are based on the variogram cloud or a binned version of it, and do not require further assumptions on the fields (Cressie 1993).For example, the package gstat relies mainly on these methods.The second often used method is the maximum likelihood (ML) method which relies on the additional important assumption of the Gaussianity of the field.Showing that the ML approach is not too far from an established approach in approximation theory, Scheuerer (2011) gives some indication why ML often works well even when the assumption of Gaussianity is violated.The ML method and its variants (geoR, Ribeiro Jr and Diggle 2001, CompRandFld, Padoan and Bevilacqua 2013) is frequently implemented for parametric inference on spatial data.For Gaussian random fields, this amounts to estimating a parametrized trend and the covariance structure.Therefore, ML methods can be readily extended to the multivariate case.
At least three different kinds of well known difficulties exist with the classical ML approach: first, the number of data is rather limited as the multivariate Gaussian density function involves the inverse of the covariance matrix, so that the calculation of the likelihood needs of order N 3 operations.There are various work-arounds, for instance, Markov random field theory allows to work directly on the inverse matrix, i.e., the precision matrix, (Rue 2001;Simpson, Lindgren, and Rue 2011).Composite likelihood approaches simplify the likelihood function by additional assumptions, cf.Lindsay (1988).A second difficulty is that the number of parameters increases rapidly with the number m of components resulting frequently in likelihood functions with lots of local maxima.As numerical optimizers such as optim are not always reliable in this case, RandomFields estimates a sequence of simpler models, starting with the model that has independent components, ending up with the full set of parameters, see Section 5 for an example.The price to be payed are rather long running times.
The third difficulty is the choice of the initial values of the numerical optimization algorithm for maximizing the likelihood.Various strategies exist.In geoR, initial starting values can be created by various methods and are then passed as an argument to the fitting procedure.
Others seem to use fixed values only (fields).Data based initial values are considered by CompRandFld using a weighted moment estimator and by gstat extracting values from the variogram cloud.Usually, the packages allow for an additional starting value given by the user.RandomFields estimates variances and scales very roughly from the data and uses fixed values for other parameters to perform several variants of the least squares fit.Afterwards it uses the best parameter set of the least squares fits as initial values for the ML estimation.Similar to the packages CompRandFld and fields, the package RandomFields recognizes sparse matrices and uses the package spam (Furrer and Sain 2010) for calculating the likelihood.
Being also mathematically based on the inverse of the covariance matrix, spatial prediction methods face similar problems as the ML approach (Chilès and Delfiner 1999).Work-arounds that share the same idea with the composite likelihood approach, such as the neighbourhood kriging (De Marsily 1984), have been suggested.Many packages exist that deal with kriging, for instance, constrainedKriging (Hofer and Papritz 2011), DiceKriging (Roustant, Ginsbourger, and Deville 2012) and geoR.The package gstat offers an enormous amount of variants for kriging, including multivariate kriging for linear models of coregionalization (Goulard and Voltz 1992).RandomFields allows for basic kriging tasks, such as simple kriging, ordinary kriging and universal kriging, on a large variety of models.
A last, also well-known difficulty appears with earth data as only few packages like sspline (Xie 2013) handle spherical data.Usually, the data are transformed in different ways to more convenient distances by the various packages.The packages gstat and fields use great circle distances.Additionally, fields also considers the data as having fully 3-dimensional coordinates, from which the usual Euclidean distances are calculated.The latter is also the case in RandomFields.While the package fields approximates the earth by a ball, RandomFields assumes an ellipsoid.

Models
Among the covariance models considered as the most important ones in geostatistics are the Whittle-Matérn model (Stein 1999;Guttorp and Gneiting 2006) support (Chilès and Delfiner 1999).Only a few packages allow for a wider spectrum of models, for instance geoR.On the multivariate side, the packages gstat and spBayes allow for some choice of multivariate models.
In the following, we describe the cross-covariance models and the cross-variogram models that can be used within RandomFields.A comprehensive list of models is given by Table 1.
Multivariate models like the latent dimension model and the linear model of coregionalization can be derived from a wide spectrum of univariate models, some of them being presented along this section.
The pieces of code given suffice to produce all the pictures as presented if the following three options R> RFoptions(seed = 0, height = 4, always_close_screen = FALSE) are set once, at the very beginning.The first option sets the seed for every subsequent simulation to 0. The second one sets the height of a single figure in the graphical display to 4 and indicates also that the aspect ratio of the figure will be calculated internally.The third one has the effect that a new window is opened for each plot.

Multivariate models
The most commonly used model for a multivariate process Z = (Z 1 , . . ., Z m ) is probably the linear model of coregionalization (Goulard and Voltz 1992), which is also offered by gstat.
The cross-covariance function of Z equals M CM where C is a diagonal matrix containing the covariance functions of the univariate processes Y j .Obviously, the vector of independent processes Y = (Y 1 , . . ., Y k ) might be replaced by any arbitrary second order process Y , so that C can be replaced by any arbitrary k × k cross-covariance function.The function RMmatrix in RandomFields allows for this generality.The parameter k must be at most 10.To consider an example, let k = 2, and Y 1 , Y 2 be independent univariate Gaussian random fields with Whittle-Matérn covariance function Here, K is a modified Bessel function, Γ the Gamma function, and ν > 0 a smoothness parameter.The values of M 1 and M 2 imply that all the variables have variance 1.This model can be defined in RandomFields as follows, see Figure 1 for a realization: R> M1 <-c(0.9,0.6) R> M2 <-c(sqrt(0.19),0.8) R> model <-RMmatrix(M = M1, RMwhittle(nu = 0.3)) + + RMmatrix(M = M2, RMwhittle(nu = 2)) R> x <-y <-seq(-10, 10, 0.3) R> simu <-RFsimulate(model, x, y) R> plot(simu) In particular, ν takes here the value 0.3 and 2, respectively.Note that, in RandomFields, all commands that define models start with 'RM' whilst functions start with 'RF'.The models in RandomFields mostly include more parameters than given here.These parameters are set to standard values if they are not given.For instance, the standard value for variance and scale is 1.Wackernagel (2003) suggests a bivariate delay effect model of the form where C 0 is an arbitrary univariate stationary covariance function and s ∈ R d .Here, the components of the corresponding random fields experience a pure shift, see Figure 2. The corresponding code is R> model <-RMdelay(RMstable(alpha = 1.9, scale = 2), s = c(4, 4)) R> simu <-RFsimulate(model, x, y) R> plot(simu, zlim = 'joint')  Here, the shift equals s = (4, 4).The underlying Gaussian random field possesses the so-called stable covariance model, also called power exponential model, with scale parameter S = 2 and shape parameter α = 1.9.The parameter zlim that is well-known from the function image allows here for further options.Both, the value 'joint' and a simple two-dimensional vector fix the same range for all images.A matrix-valued zlim is able to determine a separate range for each column of the plot.
A third basic construction principle is the latent dimension model (Apanasovich and Genton 2010).Here, additional spatial dimensions serve as a source for the multivariate components: for a random field Z 0 in R d+k and vectors z i ∈ R k .An example is given in Figure 4, where Here, the covariance function of Z 0 is a generalized Cauchy model (Gneiting and Schlather 2004), with parameter α = 1.5 and β = 3.The argument MARGIN.slicesgives the latent dimension d + k.The argument n.slices equals m if we are interested in a m-variate field.The slices are equally distributed on a grid given for the latent dimension.As this plotting technique of slices is independent of any theory on latent dimensions the names of the arguments refer to the technical aspect only.The random fields for the different variables are visually similar for all models given here.In particular, the paths of all components show the same degree of smoothness.The multivariate quasi-arithmetic mean model by Porcu, Mateu, and Christakos (2009), for some ρ red ∈ [−1, 1] and Gneiting et al. (2010) show that the right-hand side of the above equality is 0, i.e., c 12 is necessarily 0 if and only if ν 12 < 1 2 (ν 11 +ν 22 ).Hence, we may introduce ν red = 2ν 12 /(ν 11 +ν 22 ) which must be greater than or equal to 1 to have dependent components.Figure 5 gives an illustration obtained by the following code: i.e., in the above model we have For the general m-variate case, Gneiting et al. (2010) give a sufficient condition for C being a cross-covariance model, Here, RMschur defines the Schur product between two matrices, here ρ and C(h) for any fixed value of h.
Finally, Bevilacqua et al. (2013) derive multivariate models with compact support based on the Gneiting class.The bivariate model is available as RMbigneiting.

Vector-valued models
Vector fields are used to describe fluxes of liquids or gases in fluid mechanics (Patankar 1980;Chorin and Marsden 1990), as well as fields of forces, e.g., electro-magnetic fields (Stratton 2007) or, more generally, as gradients of potentials (Kellog 1953).
Stochastic vector fields appear naturally within Kolmogorov's theory of turbulent flows (Pope 2000, Section 6.2).Here, micro-fluctuations are described by an intrinsically stationary, isotropic vector field, given by the variogram γ, Here δ ij = 1 for i = j and δ ij = 0, otherwise.It is well-known that the marginal distributions in a turbulent flow are not Gaussian.
Let F = (F 1 , . . ., F d ) be a vector field on R d .The curl ∇ × F is defined as ∇ × F = ∂F 2 /∂x 1 − ∂F 1 /∂x 2 if d = 2 and as Here, the e i are the canonic unit vectors.The divergence ∇ • F is given as An important statement, not only for the theory of turbulence, is the Helmholtz decomposition: Under mild regularity assumptions, a vector field F in three dimensions can be uniquely decomposed into a curl free vector field F c , i.e., ∇ × F c = 0, and a divergence free vector field F d , i.e., Sagaut and Cambon 2008).A curl free vector field F c can either be interpreted as the flux of an incompressible liquid that is verbally driven by sinks and sources, or it describes a fluid that is compressible and liquid 'disappears' if it is compressed.The absolute value of the divergence ∇ • F gives the strength of the sinks and sources, where the sign determines whether it is a sink or a source.Divergence free vector fields describe vortices and magnetic fields.The direction of the curl ∇ × F gives the axis of the vortex and the norm of the curl its strength.
Stationary and non-stationary Gaussian random fields that are divergence free or curl free (in the mean square sense, hence pathwise) can be easily constructed from a sufficiently often differentiable univariate covariance function, C 0 say, cf.Matheron (1972Matheron ( , 1973)), for instance.For simplicity, let C 0 be translation invariant.Then, is the cross-covariance function of a stationary divergence free random field and is the cross-covariance function of a curl free random field.Here, E denotes the identity matrix and ∆ denotes the Laplace operator.Scheuerer and Schlather (2012) show that if a vector-valued random field on R 2 is divergence free or curl free or a vector-valued random field on R 3 is curl free, then their cross-covariance functions have been obtained by essentially the above construction.
The curl free models, the divergence free models, and Kolmogorov's cross-variogram have been implemented in RandomFields.The covariance model RMcurlfree defines random fields in R 2 that have four components: the first component is the (random) potential field given by C 0 , the next two components give the derived (curl free) vector field, whereas the last component gives the field of sinks and sources of the vector field, i.e., the divergence of the vector field.In fact, due to the included field of sinks and sources, whose covariance function equals ∆ 2 C 0 , the covariance function C 0 must be at least four times differentiable.Random fields simulated with RMdivfree have also four components that can be interpreted in a similar way, the last component being here the curl of the respective divergence free vector field.Figure 7 illustrates such a four-variate random field.The code is R> model <-RMcurlfree(RMmatern(nu = 5), scale = 4) R> simu <-RFsimulate(model, x, y) R> plot(simu, select.variables= list(1, 2 : 3, 4)) R> plot(model, dim = 2, xlim = c(-3, 3), main = "", cex = 2.3, col = "brown") Here, plot has the additional argument select.variables,by which the user chooses the components to be plotted and the way they are plotted.If select.variables is not given then plot plots all components separately.If select.variables is a list, then any component of the list should be a scalar or a vector and is plotted in a separate figure.If an element of the list is a scalar, an image is created.If the element is a bivariate vector a vector field is created and if the element is a trivariate vector, an image is plotted for the first component, and, additionally, a vector field is shown given by the subsequent two components.
The model RMvector allows combinations of a divergence free and a curl free spatial or spacetime random field where the ∇-operators and the Laplace operator act on the spatial components only and a ∈ [−1, 1].The model is parametrized in a way such that a = −1 corresponds to a curl free field and a = 1 to a divergence free field.The corresponding random field is bivariate and contains the vectors only.Here, the covariance function C 0 must be only twice differentiable in the spatial component since the field of sources and sinks or curls is not included.Finally, Figure 8 shows an intrinsically stationary Gaussian vector field with Kolmogorov's variogram, which has been obtained by R> x <-y <-seq(-2, 2, len = 20) R> model <-RMkolmogorov() R> simu <-RFsimulate(model, x, y, z = 1) R> plot(simu, select.variables= list(1 : 2), col = c("red")) R> plot(model, dim = 3, xlim = c(-3, 3), MARGIN = 1 : 2, cex = 2.3, + fixed.MARGIN = 1.0, main = "", col = "brown") Here, MARGIN gives the two coordinate directions that should be plotted if the dimension d is greater than 2. The argument fixed.MARGIN gives the values of the remaining coordinates, here one coordinate, that are kept fix.In fact, fixed.MARGIN = 0.0, which is the default, does not yield very illuminating results, cf.Equation 4, so that we have chosen fixedMARGIN = 1.0.

Example: weather data
Figure 7: Simulation of a curl free vector field.Top left: potential field (variable 1); top center: the corresponding curl free vector field (variables 2/3) is given by the arrows whose length can be visibly inspected by comparing them to the above arrow which represents a vector of length 2; top right: the field of sinks and sources (variable 4), i.e., the divergence of the vector field in the top center; bottom: contour plots of the covariance function for these variables 1-4.The data set considered here gives the difference between both the pressure and the temperature at 157 locations in the United States measured at a certain time instance and their predictions 48 hours ago.As in Gneiting et al. (2010) we assume that the expectation of the differences is zero, and refer to data(weather) for details.
We consider first a simplified version of the bivariate Whittle-Matérn model from Section 4.1 where ν red is set to 1 and the scale s 11 = s 12 = s 22 is the same for all entries in the cross covariance.A bivariate nugget effect with independent components is added: This model is called 'parsimonious model' in Gneiting et al. (2010).We transform the earth coordinates to Euclidean distances before performing the analysis: R> data(weather) R> Dist.mat 3 : 4])) The function RFfit recognizes all parameters that have the value NA, and estimates them through ML.Fitting the parsimonious model yields Here, the null model is the simple multivariate model and the alternative model is the parsimonious model.The RFratiotest reports the test statistic only for these two models, as the parsimonious model and the independent model do not include each other.The latter is true since the parsimonious model is more flexible with regard to the correlation structure while the independent model is more flexible with regard to scaling, cf. the arguments c and s, respectively.Following geoR, the cross validation R> RFcrossvalidate (pars, x = weather[ , 3 : 4], data = PT, full = TRUE) fits the model only to the complete data set and then the value at each of the 157 locations is predicted by kriging given the data at all the other locations.Again, the argument full=TRUE produces information about the internal sequence.The reported deviations are about the same for all three models.Hence, cross validation does not give a clear preferance to any of the three models.
In Gneiting et al. (2010), also the whole bivariate Whittle-Matérn model has been considered, cf.Section 4.1.Re-analysis of this model can be done in a similar way as for the parsimonious model, yields that the parsimonious model should be preferred.
Finally, Figure 9 shows the kriged field for the parsimonious model that interpolates pressure and temperature.To this end, the earth coordinates have been projected onto the plane.

Concluding remarks
Here, the multivariate features of RandomFields version 3.0 have been presented.Besides the extension of the tool box for analyzing and simulating random fields, a fundamental goal of version 3.0 of RandomFields has been the improvement of the accessibility of the package.To this end, the number of functions and the number of visible arguments have been reduced, at the same time as the functionality of the package has been extended.Models are now defined as functions, and even the "tilde" definition (~) of models can be used, in accordance with the formula specification of linear models.The number of man pages has been largely increased, offering a separate man page for each model.Further, several summary pages are included, and the documentation of more advanced features is shifted to secondary pages.
Currently, RandomFields offers more than 35 simple models and more than 30 functionals operating on simple models to construct more complex ones.Besides multivariate models also space-time models are included.RandomFields offers a modular construction system that allows for any combination of models, functionals operating on models, and functions, whenever the combination is mathematically sound.Validity checks of the model are performed by RandomFields.
The inference on Gaussian random fields will be extended towards mixed models with a spatial Gaussian component.The package will also advance into the analysis and the modeling of non-Gaussian fields, in particular max-stable random fields.
In Section 5 all the bivariate data have been co-located.This is certainly not the case for many data sets.Now, co-located multivariate data ease algorithms for estimation and spatial prediction.RandomFields is being developed towards arbitrary data where missing values are coded by NAs.

Figure 1 :
Figure 1: Simulation of a bivariate linear model of coregionalization.

Figure 3 :
Figure 3: Bivariate generalized delay effect model according to Equation 2: simulation (top) and contour plots for the components of the corresponding matrix-valued covariance function (bottom).

Figure 4 :
Figure 4: Simulation of the latent dimension model.
The multivariate model for three or more variables is available in RandomFields by the function RMparswm.The parsimonious model K introduced inGneiting et al. (2010) allows for more flexibility, K ij = ρ ij C ij where C is given by the equalities in Equation 3 and ρ is an arbitrary covariance matrix.This generalization is implemented by the model RMparswmX defined as R> RMparswmX <-function(nudiag, rho, var, scale, Aniso, proj) { + RMschur(M = rho, RMparswm(nudiag, var, scale, Aniso, proj)) + }

Figure 6 :
Figure 6: Simulation of an anisotropic model of coregionalization.

Figure 8 :
Figure 8: Kolmogorov's model of turbulence.The top picture shows the vector field of the first two components of a simulated three-dimensional vector field at the locations where the third coordinate is kept fixed.Their length can be visibly inspected by comparing them to the above arrow which represents a vector of length 2. The bottom pictures show the variogram values between all components at the locations where the third component is fixed to the value 1.

Figure 9 :
Figure 9: The kriged bivariate random field for the atmospheric data considered in Section 5.
Gneiting et al. (2010)l) fitted variable '' is close to or on the effective lower boundary Hence the gradient of the likelihood function might not be zero and none of the reported 'sd' values might be reliable.The main part is table below User's variables.Here, the rows deliver the estimated values and their standard deviations, respectively.Each column name indicates which parameter is estimated and consists of three parts separated by dots.Let us consider the first name 'matrix.M.1'.The left part 'matrix' refers to the model 'RMmatrix', the middle part 'M' to the parameter M of the model and the right part '1' numbers consecutively the entries of the parameter.In case the parameter is a scalar, no number is attached.The names matrix.M.1 and matrix.M.2 denote the nugget effect of the first and second component, biwm.sdenotes the scale s, biwm.nudiag.1 and biwm.nudiag.2refer to ν 11 and ν 22 , biwm.cdiag.1 and biwm.cdiag.2denotec11andc22,and biwm.rhoredrefers to ρ red of the bivariate Whittle-Matérn model in Section 4.1.As matrix.M.2 is close to zero, a note is displayed at the end and the boundaries used in optim are given, see the above table after Internal variables.Note that the values found inGneiting et al. (2010)deviate slightly, as the former and the current optimization algorithm differ.Note also that due to internal representation issues the values for the estimated standard deviation reported above differ slightly for the 'Internal variables'.As mentioned in Section 3 not only the parsimonious model is fitted, but a sequence of models.Note that the (internal) fitted variable '' is close to or on the effective lower boundary Hence the gradient of the likelihood function might not be zero and none of the reported 'sd' values might be reliable.Note that the (internal) fitted variable '' is close to or on the effective lower boundary Hence the gradient of the likelihood function might not be zero and none of the reported 'sd' values might be reliable.The first model refers to the independent model.The first two lines of the table after User's variables correspond to the first component, the last two lines to the second component.Parameters that have a dash are not estimated for the respective component.In the second model, the components are also assumed to be independent (ρ red = 0), but share the same scale parameter.For the fixed parameter ρ red , reporting the standard deviation is not appropriate, indicated by a dash.Here, the independent model, which is the simplest model concerning parameter estimation, has been fitted first.The fitted parameters serve as initial values for estimating the parameters of the second model, called simple multivariate model here.The parameters of the latter are finally used as starting values for the parsimonious model.As the parsimonious model has the smallest AIC, it should be chosen among the three models.This finding is supported by the likelihood ratio test which uses a χ 2 approximation: