# Censored pairwise likelihood-based tests for mixing coefficient of spatial max-mixture models

Max-mixture processes are defined as Z = max(aX, (1 -- a)Y) with X an asymptotic dependent (AD) process, Y an asymptotic independent (AI) process and a ∈ [0, 1]. So that, the mixing coefficient a may reveal the strength of the AD part present in the max-mixture process. In this paper we focus on two tests based on censored pairwise likelihood estimates. We compare their performance through an extensive simulation study. Monte Carlo simulation plays a fundamental tool for asymptotic variance calculations. We apply our tests to daily precipitations from the East of Australia. Drawbacks and possible developments are discussed.

## Authors

• 2 publications
• 10 publications
• 1 publication
• ### A Model-Free Selection Criterion For The Mixing Coefficient Of Spatial Max-Mixture Models

One of the main concerns in extreme value theory is to quantify the depe...
01/03/2018 ∙ by Abul-Fattah Abu-Awwad, et al. ∙ 0

• ### A Bayesian semi-parametric hybrid model for spatial extremes with unknown dependence structure

The max-stable process is an asymptotically justified model for spatial ...
03/24/2020 ∙ by Yuan Tian, et al. ∙ 0

• ### A semi-parametric estimation for max-mixture spatial processes

We proposed a semi-parametric estimation procedure in order to estimate ...
10/23/2017 ∙ by M. Ahmed, et al. ∙ 0

• ### Estimating a mixing distribution on the sphere using predictive recursion

Mixture models are commonly used when data show signs of heterogeneity a...
10/20/2020 ∙ by Vaidehi Dixit, et al. ∙ 0

• ### Maximum likelihood estimation for a general mixture of Markov jump processes

We estimate a general mixture of Markov jump processes. The key novel fe...
03/03/2021 ∙ by Halina Frydman, et al. ∙ 0

• ### Max-sum tests for cross-sectional dependence of high-demensional panel data

We consider a testing problem for cross-sectional dependence for high-di...
07/08/2020 ∙ by Long Feng, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The rise of risky environmental events leads to renewed interest in the statistical modelling of extremes, for example modelling extreme precipitation is pivotal in flood protection. In the last decade, max-stable (MS) models have arised as a common tool for modeling spatial extremes, since they extend the gereralized extreme value (GEV) distribution to the spatial setting, providing a consistent multivariate distributions for maxima in arbitrary dimensions.

Max-stable (MS) processes for spatial data were first constructed using the spectral representation of [13]. Subsequent developments have been done on the construction of MS process models, [27], [26], [17], and [9].

Despite many attractive properties of MS models, these processes are restricted since only AD or exact indpendence can be modeled. This drawback is constraining when modeling the tail behavior of the multivariate distribution of the data, since it is difficult to asses in practice whether a data set should be modeling using asymptotic dependent (AD) or asymptotic independent (AI) models [28] and [10]. In [32], a flexible class of models which may take into account AD and AI is proposed. Theses models, called max-mixture (MM), are a mixture of a MS process and an AI process: , with an AD process, an (AI) process and . So that, represents the proportion of AD in the process and our purpose is to propose statistical tests on the value of .

To the best of our knowledge, only a few papers address testing AI issue for spatial extreme fields. E.g., in [1]

a statistical test for AI of a bivariate maxima vector is proposed and a generalization to spatial context is done. It is based on the

- madogram [6].Standard likelihood inference on the parameters of MS models is not possible in general, since the full likelihood is not easily computable for MS vector in dimension greater than . Composite likelihood are usefull when the fully specified likelihood is computationally cumbersome, or when a fully specified model is out of reach [21] and [29]. Maximum pairwise likelihood estimation for MS models weren first suggested by [24] and is now widely used. In particular, in [32] and [2] the parameter inference for MM processes have been driven.

In this paper, we propose a test on the mixing value of a MM process. This is achieved using pairwise likelihood statistics. The paper organised as follows. Section 2 reviews the theory of spatial extreme processes MS and MM. The censored pairwise likelihood approach is presented for the statistical inference in Section 3, our testing proposal approach and the main properties are detailed in Section 4. In Section 5 we show, by means of a series of simulation studies the performance of our proposed tests. In section 6 we illustrate our testing approach by the analysis of daily precipitation from the East of Australia. Concluding remarks and some perspectives are addressed in Section 7.

## 2 Spatial extremes modeling

Throughout our work, , (generally, ) is a spatial process will be assumed to be stationary and isotropic.

### 2.1 Max-stable processes

Let , are independent replicates of a stochastic process . Then is a MS process if and only if there exist a sequence of continuous functions and such that the rescaled process of maxima, , converges in distribution to (see [12] for more details). By this definition, MS processes offer a natural choice for modeling spatial extremes.

The univariate extreme value theory, implies that the marginal distributions of are Generalized Extreme value (GEV) distributed, and without loss of generality the margins can transformed to a simple MS process called standard Fréchet distribution, .

Following [13] and [26], a simple MS process has the following representation

 X(s)=maxk≥1Qk(s)/Pk,s∈X. (2.1)

where are independent replicates of a non-negative stochastic process with unit mean at each , and are the points of a unit rate Poisson process .

The joint distribution function of the process

at locations is given by

 P(X(s1)≤z1\/,…\/,X(sK)≤zK)=Gs1\/,…\/,sK(z1,...,zK)=exp{−Vs1\/,…\/,sK(z1,...,zK)} (2.2)

where is called the exponent measure. It summarises the structure of extremal dependence and satisfies the property of homogenity of order and . It has to be noted that

 P{X(s1)≤z,..,X(sK)≤z}=exp{−Vs1\/,…\/,sK(1,...,1)/z}=exp{−1/z}θs1\/,…\/,sK\/,

with . The coefficient is known as the extremal coefficient. It can be seen as a summary of extremal dependence with two boundary values, complete dependence , and complete independence,

. In the bivariate case, the AI and AD between a pair of random variables

and , with marginal distributions and , may be identified by

 χ=limu→1−P(F1(Z1)>u|F2(Z2)>u))\/. (2.3)

The cases and represent AI and AD, respectively, [16]. This coefficient is related to the pairwise extremal coefficient through the relation .

Since both dependence functions and are useless for AI processes, [5] proposed a new dependence coefficient which measures the strength of dependence for AI processes:

 ¯χ=limu→1−¯χ(u)=limu→1−2logP(F(Z(s))>u)logP(F(Z(s))>u,F(Z(s+h))>u)−1 (2.4)

AD (respectively AI) is achieved if and if (resp. ).

Different choices for the process in (2.1) lead to some useful MS models, commonly used choices are the Guassian extreme value process [27], the extremal Gaussian process [26], the Brown-Resnick process [17], and the extremal t process [22]. Below, we list two specific examples of MS process models.

The storm profile model [27], is defined by taking , where is a Guassian density with covariance matrix , and

is a homogenous Poisson process. The bivariate marginal probability distribution of the Smith model

has the form

 P(X(s)≤z1\/,X(s+h)≤z2)=Gh(z1,z2)=exp{−[1z1Φ(γ(h)2+1γ(h)logz2z1)+ 1z2Φ(γ(h)2+1γ(h)logz1z2)]}

where , in which is the seperation vector between the two locations,, and , and

is the standard normal distribution. The pairwise extremal coefficient is

.

The Truncated Extremal Gaussian (TEG) model is originally due to [26], it is defined using

 (2.5)

where are independent replicates of a stationary Gaussian process with zero mean, unit variance and correlation function . is the indicator function of a compact random set , are indepndent replicates of and are points of Poisson process with a unit rate on . The constant is chosen to satisfy the constraint .

The bivariate marginal probability distribution of TEG model in the stationary case has the form

 P(X(s)≤z1\/,X(s+h)≤z2)=Gh(z1,z2)=exp{−(1z1+1z2)⋅ (2.6) [1−α(h)2(1−√1−2(ρ(h)+1)z1z2(z1+z2)2)]}

where is the spatial lag, and if is a disk of fixed radius . The pairwise extremal coefficient . TEG model were fitted by [9] to extreme temperature data.

The stochastic process defined in equation(2.1) has the bivariate density function

 g(z1,z2)={∂∂z1VX(z1,z2)∂∂z2VX(z1,z2)−∂2∂z1z2VX(z1,z2)}exp{−VX(z1,z2)}

where is the exponent measure of the MS process .

### 2.2 Hybrid models of spatial extremal dependence

Although MS models seems to be suitable for modeling extremely high threshold exceedances, AI models may show a better fit at finite thresholds. Since it may be difficult or impossible in practice to decide whether a dataset should be modeled using AD or AI, [32] have introduced an hybrid spatial dependence model, which is able to capture both AD and AI.
Let , be a stationary simple MS process, and , be a stationary AI process with unit Fréchet margins (see below for the construction of such a process). Then for a mixture proportion , a spatial max-mixture (MM) process is constructed:

 Z(s)=max{ aX(s),(1− a)Y(s)} (2.7)

The bivariate distribution function for a pair of sites is straightforwardly obtained by the independence between the processes and

 Gh(z1,z2)=GXh(z1 a,z2 a)GYh(z11− a,z21− a) (2.8)

where (resp. ) is the bivariate distribution function of , with space lag (resp. of ).
AI processes with unit Fréchet marginal distributions can easily be constructed (see [32]). Consider
, where is a Gaussian process. Then, is an AI process with unit Fréchet marginal distributions. Another class of AI processes called inverted max-stable (IMS) processes are defined using a simple MS process , let

 Y(s)=−1/log{1−exp[−X(s)−1]} (2.9)

With this construction, any MS process may be transformed to provide an AI counterpart. Bivariate distribution function and density of the margins of are respectively

 GYh(z1,z2)=−1+exp(−z−11)+exp(−z−12)+exp{−VYh[ω(z1),ω(z2)]} (2.10)
 gYh(z1,z2) ={∂∂z1VYh[ω(z1),ω(z2)]∂∂z2VYh[ω(z1),ω(z2)]−∂2∂z1z2VYh[ω(z1),ω(z2)]} ×exp{−VYh[ω(z1),ω(z2)]}ω2(z1)ω2(z2)exp{−(z−11+z−12)}z21z21{1−exp(−z−11}{1−exp(−z−11}

where is the exponent measure of the bivariate distribution , and .
Thus, in the case where is a MS process and is a IMS process, the distribution function in (2.8) has the form

 exp{−aVX(z1,z2)}{−1+exp[−(1−a)/z1]+exp[−(1−a)/z2]+exp[−VY[ω(1−a/z1),ω(1−a/z2)]} (2.11)

[2] analyzed daily rainfall data in the east of Australia with a class of different models (MS, AI, and MM), and showed that MM models has the merit to overcome the limits of MS models in which only AD or exact independence can be modeled.

## 3 Inference for MM processes: censored pairwise likelihood approach

In order to propose a testing procedure on on the mixing coefficient a of MM processes defined by equation (2.7), we shall use the composite likelihood.

The composite likelihood technique [21] is a general method of inference for dealing with large datasets and/or miss-specified models. A composite likelihood consists of a combination of valid likelihood objects usually related to small subsets of data and defined as

 LC=K∏k=1g(z∈Bk;φ)ωk (3.1)

where is a paremetric statistical model, is a set of marginal or conditional events, is a set of suitable weights, if the weights are all equal they may be ignored, non-equal weights may be used to improve the statistical performance in certain cases. The associated composite log-likelihood is .

In the spatial setting, the definition of a pairwise log-likelihood is derived from (3.1) by taking as the set of bivariate subvectors of taken over all disitinct sites pairs and . Thus the weighted pairwise log-likelihood is given by

 pℓ(ϑ;z)=M∑k=1K−1∑j=1K∑j′=i+1ωjj′logL(zk(sj),zk(sj′);ϑ)) (3.2)

where are the data available on the whole region, is the th observation of the th site, and is the likelihood function based on observations at locations , and , are non negative weights specifying the contributions of each pair. A simple weighting choice is to let , where is the pariwise distance, and is a specified value.

Inference using pairwise likelihood methods is computationaly expensive, since with sites there are pairs to include. This methodology has been used by [24], [9] and [11] for infernence on MS processes.

Different inference approaches based on a censored threshold-based log-pairwise likelihood have been used by several researchers [20], [32], [2] and [15]. Where the censored pairwise contributions take the forms

 Lu(zk(sj),zj(sj′);ϑ)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩G(u,u;ϑ),if max (zk(sj),zk(sj′))≤u∂zk(sj)G(zk(sj),u;ϑ),if zk(sj)>u,zk(sj′)≤u∂zk(sj′)G(u,zk(sj′);ϑ),if zk(sj)≤u,zk(sj′)>ug(zk(sj),zk(sj′);ϑ),if max (zk(sj),zk(sj′))≥u (3.3)
 Lu(zk(sj),zj(sj′);ϑ)={G(u,u;ϑ),if max (zk(sj),zk(sj′))≤ug(zk(sj),zk(sj′);ϑ),if max (zk(sj),zk(sj′))≥u (3.4)

where is a high threshold, is given in equation (2.8) and is the bivariate density function, i.e. .

## 4 Pairwise likelihood statistics for testing H0:a=a0 versus H1:a≠a0

We assumed the parameters of a MM model is partitioned as , where is a one-dimensional parameter of interest that denotes the mixing coefficient for a MM model and is a nuisance parameter, .

Our purpose is to test the hypothesis versus , for some specified value [0,1]. Let denotes the unrestricted maximum pairwise likelihood estimator and , denotes the constrained maximum pairwise likelihood estimator of for a fixed . The pairwise maximum likelihood estimator is asymptotically normally distributed:

 (4.1)

where is the -dimensional normal distribution with mean and variance , and denotes the Godambe information matrix:

 G−1(ϑ)=H−1(ϑ)J(ϑ)H−1(ϑ)\/,

where is called the sensitivity matrix and is called the variability matrix. For more details see [18], [21], and [30].

In this paper we propose the following two statistics exploiting the pairwise maximum likelihood as an inferential tool. Our objective to facilitate the modeling of the spatial data by a random field with appropriate extremal behaviour. The

-test statistic which is

straightforwardly derived from thecentral limit theorem for maximum composite likelihood estimators.

 Z=^a−a√Gaa(^ϑpℓ)D−→N{0,1} (4.2)

where denotes a submatrix of the inverse of pertaining to . While the pairwise likelihood ratio statistic (

) with a nonstandard asymptotic chi-squared distribution is given by

 LR=2{pℓ(^ϑpℓ)−pℓ(^ϑ∗pℓ}D−→λχ21 (4.3)

where , and are respectively submatrices of the inverse of and pertaining to and

is a chi square random variable with one degree of freedom.

The asymptotic distribution of the statistic has been studied in [18] in a more general context (specifically when the dimension of the parameter of interest may be greater than ). Different kind of adjustments were proposed to recover an asymptotic chi square distribution in [4], [14], [19], [23] and [25].

Standard errors and critical values for the tests require the estimation of the Godambe matrix and its components, and since analytical expressions for and are difficult to obtain in mostly realistic applications. They are usually estimated by means of a Monte Carlo simulations:

 ^HS(ϑ)=−M−1M∑k=1∇2pℓ(^ϑpℓ;zk)
 ^JS(ϑ)=M−1M∑k=1∇pℓ(^ϑpℓ;zk)∇Tpℓ(^ϑpℓ;zk)\/,

where , , denote the th datasets simulated from the fitted model. The results obtained by [3] for testing the paremeters of equicorrelated multivariate normal model showed that the coverage of the statistics based on Monte Carlo simulation are almost indentical to those of statistics based on analytically computed quantities.

Testing at boundary points (AD) or (AI) is non standard since there are additional nuisance parameters which are present only under the alternative [7] and [8], we apply our test at some points close to the boundaries, i.e. or .

## 5 Simulation study

We have performed several of simulation studies in order to investigate the performance of our testing procedures. We simulated from a MM model (2.7) in which is a TEG process (2.6) with a disk of fixed radius . The AI process is an inverse TEG process with a disk of fixed radius . For simplicity, we assume that the correlation functions of these two processes are exponential, with range parameters respectively. Our purpose is to test versus , varies from 0.01 to 0.99 by steps of 0.01.

The censored pairwise likelihood approach (3.4) is used for estimation, where the threshold

is taken corresponding to the 0.9 empirical quantile at each site, and equal weights are

considered. To reduce computational burden the pairwise likelihood function has been coded in C; the optimization has been parallelized on 20 cores using the R library parallel and carried out using the Nelder–Mead optimization routine in R.

We have done replications of independent copies of the considered MM process on locations uniformaly chosen in the the square . The Boxplots of the estimated parameters on the samples are displayed in Figure 1. The parameters used are and different mixing coefficient and are considered. Generally, the parameters are well estimated.

In order to obtain accurate estimates of and in the Monte Carlo procedure, we perform an exploratory study with 200 simulation replicates based on MM model described above with parameters . For each replication, we randomly generate locations uniformaly in the square [0, 2] [0, 2] with 1000 independent observations at each sampled location. Then, we simulate data from the fitted model with independent simulations at the sampled locations. In Figure 2, we present the boxplots for and . The results give a justification to use as a compromise between accuracy and computation time.

Figure 3 and Tables 1, 2 and 3 report a summary of comparison results between the two statistics and in terms of empirical probabilities for concluding in testing hypothesis against based on 150 simulation replicates from 1000 independent copies simulated at 50 randomly and uniformaly sampled locations in the square from three MM models with parameters , according to different values of the mixing parameter . Decisions obtained at three significance levels .

Despite the poor performance which may be expected at the region around to the true mixing coefficient . The results in terms of empirical probabilities for concluding of the two statistic show a reasonable performance. The probability of making a correct decision becomes higher as we become very close or move far from the model true value. We also remark that the performances of the two tests are very similar except for , where the test presents an unexpected over sensitivity.

We have also explored the AD and AI cases (i.e and respectively).The two tests were performed for all values in the set . We have computed the probability of concluding while the true parameter is (AD case) or (AI case). For this purpose a MM model were fitted for each 150 simulation replicates from a MS TEG and inverted TEG processes with parameters , , repectively, and a moderately sized data with = 50 sites and =1000 independent observations. The spatial domains respectively are (AD case) and (AI case). The results are presented in Figure 4 and Tables 4 and 5.

Despite the interesting results, that show the increase of the empirical probability of concluding , both tests fail to identify precisely AD or AI. They could be used as indication of AD or AI but are not sensitive.

## 6 Rainfall Data Example

The data analysed in this section are daily rainfall amounts in (millimetres) over the years 1972-2016 occurring during April-September at 38 sites in the East of Australia whose locations are shown in Fig. 5. The altitude of the sites varying from 4 to 552 meters above mean sea level. The sites are separated by distances from km to km. This data set is freely avaliable on http://www.bom.gov.au/climate/data/.

Following the approach of [2], we have done a graphical exploration using the coefficients and , in order to evaluate a possible anisotropy in the data. Figure 6 is based on the empirical estimates of the functions and in different directional sectors , , , and , where represents the northing direction. We can conclude that there is no evidence for anisotropy and we shall consider that the data comes from an isotropic process.

A way to apply our testing approach is based on dividing the daily rainfall dataset from the 38 sites into two groups A and B (see Figure 5). We shall consider five models which belong to the three classes: MM, MS, and IMS. These models are first fitted on data from group A. The composite likelihood information criterion (CLIC) [31], defined as CLIC= is used to choose the best fitted model. Lower values of CLIC indicate a better fit. Then we apply our test with the best fitted model for group A. Let be the mixing parameter of the best fitted model for group A (we found ). Our test will be done on the group B data, and will be vs .

The fitted models are:

: a MM model where is a TEG process with an exponential correlation function , . is a disk of fixed and unknown radius , and is an inverted TEG process with exponential correlation function , , and is a disk with fixed and unknown radius .

: a MM model where is a TEG process as in . is an isotropic inverted Smith process where is a diagonal matrix () with , i.e., .

: a MS TEG process described as in .

: a MS isotropic Smith process where is a diagonal matrix () with , i.e., .

: the inverted Smith process described as in .

The considered models have unit Fréchet marginal distributions. We fit a GEV distribution on each site and then transform the marginal laws to unit Fréchet using

 ^F(z)=(1+^ξz−^μ^σ)1/^ξ\/,

where are the estimated parameters of the GEV distribution. The censored pairwise likelihood approach (3.4) is used in order to estimate the parameters. The threshold is and equal weights are used. The matrices and and the related quantities CLIC and standard errors are obtained by Monte Carlo procedure through simulating data with independent draws at the sampled 19 sites from the fitted model.

Our results are summarised in Table 6. The best-fitting model for group A, as judged by CLIC, is the hybrid dependence model . The mixing coefficient for the other hybrid model is very close to one, which indicates that there is no mixture between the max-stable process and the asymptotically independent one, so the asymptotic independence components are not identifiable. This fact affects the values of the estimates. Moreover, model reduces to model .

For data from group B, we have considered the two statistics and described in this paper, to test if the hybrid model can be used to make inference for this group, i.e., testing versus . We obtain the calculated values for and , (value), , , (value).

This leads us to retain and thus that there is no differences in the mixing parameter between the two groups. We have also performed an independent two-samples -test: let (resp. be the mixing parameter for group A (resp. group B), consider . The statistic

 ZC=^aA−^aBSE^aA−^aB

where stands for the estimated standard error. The calculated value of is with -value , leading to retain and conclude that there is no significant difference between the two mixing coefficient.
Nevertheless, these conclusions are subordinated to the assumption that both groups A and B have the same underlying model M2.

## 7 Conclusion

In this paper we have considered hypothesis testing for the mixing coefficient of a MM models proposed in [32] using two statistics the and the when a censored pairwise likelihood is employed for inferential purposes.

The two statistics has emerged as an efficient tools for testing hypothesis on the mixing coefficient with better performance achieved by , but with the drawback of a nonstandard asymptotic distribution at the boundaries, since the number of nuisance parameters is different between the two hypothesis and it also requires heavier computations. Our procedure seems to be a performant validation tool.

One other drawback of our work is that the proposed tests model-dependent. In the future, we plan to propose a free-model test, using the F-madogram.

Acknowledgments. We are grateful to Jean-Noël Bacro, Carlo Gaetan and Gwladys Toulemonde for giving their estimations C codes which we used as a base for computing the statistics of our tests. This work was supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program ”Investissements d’Avenir” (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). It was also supported by the CERISE LEFE-INSU projcet.

## References

• [1] Jean-Noël Bacro, Liliane Bel, and Christian Lantuéjoul. Testing the independence of maxima: from bivariate vectors to spatial extreme fields. Extremes, 13(2):155–175, 2010.
• [2] Jean-Noel Bacro, Carlo Gaetan, and Gwladys Toulemonde. A flexible dependence model for spatial extremes. Journal of Statistical Planning and Inference, 172:36–52, 2016.
• [3] Manuela Cattelan and Nicola Sartori. Empirical and simulated adjustments of composite likelihood ratio statistics. Journal of Statistical Computation and Simulation, 86(5):1056–1067, 2016.
• [4] Richard E Chandler and Steven Bate. Inference for clustered data using the independence loglikelihood. Biometrika, pages 167–183, 2007.
• [5] Stuart Coles, Janet Heffernan, and Jonathan Tawn. Dependence measures for extreme value analyses. Extremes, 2(4):339–365, 1999.
• [6] Dan Cooley, Philippe Naveau, and Paul Poncet. Variograms for spatial max-stable random fields. Dependence in probability and statistics, pages 373–390, 2006.
• [7] Robert B Davies. Hypothesis testing when a nuisance parameter is present only under the alternative. Biometrika, 64(2):247–254, 1977.
• [8] Robert B Davies. Hypothesis testing when a nuisance parameter is present only under the alternatives. Biometrika, pages 33–43, 1987.
• [9] AC Davison and MM Gholamrezaee. Geostatistics of extremes. In Proc. R. Soc. A, volume 468, pages 581–608, 2012.
• [10] Anthony C Davison, Raphaël Huser, and Emeric Thibaud. Geostatistics of dependent and asymptotically independent extremes. Mathematical Geosciences, 45(5):511–529, 2013.
• [11] Anthony C Davison, Simone A Padoan, and Mathieu Ribatet. Statistical modeling of spatial extremes. Statistical science, pages 161–186, 2012.
• [12] L. De Haan and T. T. Pereira. Spatial extremes: Models for the stationary case. The annals of statistics, pages 146–168, 2006.
• [13] Laurens De Haan. A spectral representation for max-stable processes. The annals of probability, pages 1194–1204, 1984.
• [14] Helena Geys, Geert Molenberghs, and Louise M Ryan. Pseudolikelihood modeling of multivariate outcomes in developmental toxicology. Journal of the American Statistical Association, 94(447):734–745, 1999.
• [15] Raphaël Huser and AC Davison. Space–time mo