Marie Chavent

Benoît Liquet

Jérôme Saracco

A SEMIPARAMETRIC APPROACH FOR A MULTIVARIATE SAMPLE SELECTION MODEL

Keywords: Sliced Inverse Regression (SIR), Multivariate SIR, Canonical Analysis, Semiparametric Regression Models, Eigen-decomposition

Most of the common estimation methods for sample selection models rely heavily on parametric and normality assumptions. We consider in this paper a multivariate semiparametric sample selection model and develop a geometric approach to the estimation of the slope vectors in the outcome equation and in the selection equation. Contrary to most existing methods, we deal symmetrically with both slope vectors. Moreover, the estimation method is link-free and distributionfree. It works in two main steps: a multivariate sliced inverse regression step, and a canonical analysis step. We establish √ n-consistency and asymptotic normality of the estimates. We describe how to estimate the observation and selection link functions. The theory is illustrated with a simulation study.

Introduction

Sample selection models (SSM) are described by two equations. A selection equation specifies the state "observed / non-observed (missing)" of the dependent variable y as a function of explanatory variables x. An outcome equation specifies the value of the dependent variable y as another function of explanatory variables

x. Numerous papers dealing with univariate SSM have been published. The adjective "univariate" refers to y ∈ R. In this paper, we focus on multivariate SSM, that is, when y ∈ R q , q > 1.

Let us first briefly review of univariate SSM. [START_REF] Heckman | Sample selection bias as a specification error[END_REF] introduced what is now regarded as the prototype selection model. [START_REF] Amemiya | Advanced econometrics[END_REF] refers to this model as the type II Tobit model:

(E1) : y * 1 = θ 1 + x ′ β 1 + ε 1 (E2) : y * 2 = θ 2 + x ′ β 2 + ε 2 (E3) : y 2 = I[y * 2 > 0] (E4) : y 1 = y * 1 y 2 (E5) : (ε 1 , ε 2 ) ′ |x ∼ N (0, Γ), Γ = σ 2 1 σ 12 σ 12 σ 2 2
where the notation I designates the indicator function. The observed variables are y 1 ∈ R, y 2 ∈ {0, 1} and x ∈ R p . Note that in this model, the explanatory variable x does not include the y variable, contrary to [START_REF] Maddala | Limited-dependent and qualitative variables in Econometrics[END_REF] who considered a more general simultaneous equation modelling framework where the outcome y can appear on both right and left hand sides of the equations (E1) and (E2). Note also that, in equation ( E4), missing values are denoted by zeoro, leading to possible confusions with zero as an actual observed value for y 1 . Equation (E3) is the selection equation and equation (E2) is the potential outcome equation. The maximum likelihood method is generally used to estimate such models. The score function is highly non-linear. The convergence of the algorithm heavily depends on the choice of good initial values, and the asymptotic properties of the estimate are very sensitive to the model specification. This has been discussed by [START_REF] Goldberger | Abnormal selection bias[END_REF] among others. Alternative methods have been designed. [START_REF] Heckman | Sample selection bias as a specification error[END_REF] proposed a two-step method estimating first the selection equation, and then using the result to estimate the outcome equation in a second stage. Many authors have considered parametric estimation methods.

For a survey of these aspects, see [START_REF] Amemiya | Advanced econometrics[END_REF], [START_REF] Maddala | Limited-dependent and qualitative variables in Econometrics[END_REF][START_REF] Maddala | Estimation of Limited-dependent Variable Models under Rational Expectations[END_REF] or [START_REF] Blundell | Simultaneous microeconometric models with censored or qualitative dependent variables[END_REF].

Semiparametric estimation methods have been developed to bypass the sensitivity to specification assumptions. They handle more general models, especially for error specification. [START_REF] Melenberg | Semi-parametric estimation of the sample selection model[END_REF] give an overview of the semiparametric estimation methods for SSM. Most semiparametric estimation techniques of SSM also proceed in two stages. The first gives a consistent estimate of the slope of the selection equation. The second stage works with the non-missing y only, (i) building a biased estimate of the slope of the outcome equation, and (ii) correcting for this bias with the help of the slope estimated in the first step. [START_REF] Duan | Distribution-free and link-free estimation method for the sample selection model[END_REF], [START_REF] Newey | Two-step series estimation of sample selection models[END_REF], Ahn and Powell (1992), [START_REF] Lee | Semiparametric two-stage estimation of sample selection models subject to Tobit-type selection rules[END_REF] follow such a scheme.

In this paper, we examine multivariate sample selection models (MSSM) which are a generalization of the type II Tobit model when the dependent variable y is a vector of R q . This kind of model can also be seen as a generalization of classical multivariate Tobit model defined by: y = max(y * , 0) where y * = Cx + ε, ε ∼ N (0, Γ) and C is a q × p matrix of coefficients (see for instance [START_REF] Eiswerth | Examining Post-Wildfire Reseeding on Arid Rangeland: A Multivariate Tobit Modelling Approach[END_REF], for a brief presentation and an ecological application of this model).

We focus on a semiparametric MSSM by introducing unknown link functions in the selection and outcome equations in order to get a more flexible model.

Moreover, we do not assume that the distribution of the error term is a multivariate normal distribution. Like [START_REF] Duan | Distribution-free and link-free estimation method for the sample selection model[END_REF] in the univariate case, we propose a link-free and distribution-free estimation method. Contrary to most existing methods, we deal symmetrically with both slopes (of the selection and outcome equations).

In Section 2, we give a description of the semiparametric MSSM. We show in Section 3 the geometric approach to the estimation of the slopes of the outcome and selection equations from a population point of view and we give the corresponding sample version in order to obtain the slope estimators. The estimation method works in two steps (which have nothing to do with the two classical stages of the approaches mentioned above). The first one performs a multivariate sliced inverse regression (MSIR) analysis. The second step converts the MSIR indices to estimators of the slopes by means of two canonical analyses. The corresponding numerical algorithm is fast (since the method is based on only a few matrix calculus and eigen-decompositions, without need for any time-consuming iterative computations) and does not require starting values. Asymptotic properties of the slope estimators are derived in Section 4. Simulation results are reported in Section 5. Finally, concluding remarks are given in Section 6.

A semiparametric multivariate sample selection model

We consider the following semiparametric multivariate sample selection model:

for j = 1, . . . , q, y (j) =    g (j) 1 x′ 1 γ1 , ε (j) 1 if g (j) 2 x′ 2 γ2 , ε (j) 2 > 0 M V otherwise. (2.1)
where:

-The symbol M V symbolically indicates a missing (non observed) value for y (j) in order to avoid any confusion with zero as an observed value.

-The dependent variable y = (y (1) , . . . , y (q) ) ∈ R q (when each y (j) is observed) is a q-dimensional random vector. In the following, we will see that there is no need to require all values for the y (j) 's to be real.

-The functions g

(j)
1 and g

(j)
2 are unknown link functions. For the j-th component y (j) of y, g

(j)
1 is called the observation link function and g (j)

2 the selection link function.

-The variables x1 ∈ R p 1 and x2 ∈ R p 2 are subvectors of a random vector x ∈ R p , assumed to have an elliptically symmetric distribution with parameters µ = E(x)

and Var(x) = Σ. Let A k , k = 1, 2 be a p×p k matrix which selects the components of xk in x, that is: xk = A ′ k x. This matrix has exactly one "1" in each column and at most one "1" in each row, and the other elements are "0". From the definition of A k , this matrix is a full column rank matrix such that A ′ k A k = I p k . These matrices A 1 and A 2 are assumed to be known a priori. They are not chosen arbitrarily by the user, they need to be assumed based on existing theory on the exclusion of specific variables. It follows that x1 and x2 are elliptically distributed with parameters

µ k = E(x k ) = A ′ k µ, k = 1, 2 and Σ k = Var(x k ) = A ′ k ΣA k , k = 1, 2. -Let ε (j) = (ε (j) 1 , ε (j)
2 ) ′ . Let us also define ε = ε (1)′ , . . . , ε (q)′ ′ . The error term ε is a random vector independent of x with an unknown distribution.

-The parameters γ1 and γ2 are the p 1 × 1 and p 2 × 1 real unknown slope parameters. Let us also introduce γ k = A k γk ∈ R p , k = 1, 2, in order to expand γk to a p × 1 vector with zeros corresponding to the non-selected components.

Under the generality of the unknown link functions in this model, the intercepts, the vector lengths and vector signs of γ1 and γ2 are not identifiable.

Without additional assumptions, only the directions of the observation and selection slope vectors are identifiable. Then, our main purpose is to estimate the directions of the vectors γ1 and γ2 . The nonparametric estimation of g (j)

1 and g (j) 2 will also be discussed.

We will consider model (2.1) as a particular case of a more general multivariate two-index semiparametric regression model of the form

y = f (x ′ γ 1 , x ′ γ 2 , ε). (2.2)
Model (2.2) was introduced by Li (1991) when y ∈ R. Li (1991) introduced the sliced inverse regression in order to estimate the subspace of R p , spanned by the γ k 's, which is called the e.d.r. (effective dimension reduction) space. In model (2.2), since the link function f is assumed to be arbitrary and unknown, the γ k 's are not individually identifiable, while the e.d.r. space is identifiable.

Some extensions of the SIR approach to multivariate y have been studied by [START_REF] Aragon | A Gauss implementation of multivariate sliced inverse regression[END_REF], [START_REF] Li | Dimension reduction for multivariate response data[END_REF], [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF], and [START_REF] Barreda | Some extensions of multivariate sliced inverse regression[END_REF].

It is interesting to note that SIR and Pooled Marginal SIR (a multivariate SIR approach which will be used in the next section) do not require a metric structure for the outcome variable(s). Thus, MV values for the y (j) 's are easily managed.

In our context, we have to take into account extra information about the e.d.r. space, namely, structural zeros in the slopes γ 1 and γ 2 , with a link function f depending on the unknown functions g (j)

1 and g (j)

2 for j = 1, . . . , q. We now exhibit in Theorem 1 a geometrical property of this model on which the proposed approach is based. Let us define the linear subspace E = Span(γ 1 , γ 2 ) of R p . Without additional conditions, we have dim(E) ≤ 2. If γ 1 and γ 2 are linearly independent, then dim(E) = 2, and {γ 1 , γ 2 } is a basis of the e.d.r. space. In order to ensure that we are working on a two-index model (that is dim(E) = 2), let us assign the following identifiability conditions:

(i) Each vector xk , k = 1, 2, has at least an x-component not present in the other xk , k = 2, 1; such a component could be considered k-specific. (ii) At least one component of γ k among the k-specific component is non null, k = 1, 2.
Note that these identifiability conditions are stronger than the usual identifiability condition, which is that x2 contains an x-component that is not in x1 . The underlying reason for the stronger condition is that the proposed method deals symmetrically with the selection and outcome slope vectors. Knowingly, we do not make use of an important piece of information, namely that the selection probabilities depend only on one of the two index variables, x′ 2 γ2 . We now bring these conditions into a geometrical perspective. Let us consider the linear subspace

E k = Span(A k ) of R p .
Theorem 1. Under the assumptions of model (2.1) and the identifiability conditions, we have: for k = 1, 2,

E ∩ E k = Span(γ k ).
Proof. From the definition of A k , we have dim(E k ) = p k . The identifiability conditions give: (i)

E 1 ⊂ E 2 and E 2 ⊂ E 1 , and (ii) E ∩ E 1 = E and E ∩ E 2 = E.
Let us study more closely the linear subspace

E ∩ E k . Since dim(E) = 2, we have dim(E ∩ E k ) ≤ 2. From the definition of E and E k , γ k ∈ E ∩ E k and then dim(E ∩ E k ) ≥ 1. From the identifiability conditions, we get, for k ⋆ = k, γ k ⋆ ∈ E and γ k ⋆ ∈ E k , thus γ k ⋆ ∈ E ∩ E j and dim(E ∩ E k ) < 2. Finally, dim(E ∩ E k ) = 1 and E ∩ E k ⊂ R p is spanned by γ k .
We specify in the next section how to determine a basis of E and to deduce a basis E ∩ E j from a population point of view. Then we describe how to estimate the directions of γ 1 and γ 2 .

Remark 1. The full model defined in (2.1) can be interpreted as an item nonresponse model, that is the response status for each outcome measure (or survey item) is governed by a specific selection equation. We can also introduce a simplified model in terms of the type of missing data encountered. In the following simplified version of the model, we assume that the same selection equation is used for all outcomes: each selection link function g (j) 2 (.) is equal to the same link function g 2 (.). With unique error term ε 2 , the model can be written in this simple way:

y = g 1 (x ′ 1 γ1 , ε 1 ) if g 2 (x ′ 2 γ2 , ε 2 ) > 0 M V otherwise.
where the observation link function g 1 (.) takes its values (when they are observed) in R q and the error term ε 1 is a q-dimensional random vector. This model can be interpreted as a case non-response model, when the response status for multiple outcomes is clustered at the individual level: an individual either responds to all outcome measures (case response) or does not respond to any outcome measure (case non response).

Remark 2. The proposed approach can cope with the generalized two-limit selection model of the form: for j = 1, . . . , q,

y (j) =      L * (j) 1 if g (j) 2 (x ′ 2 γ2 , ε 2 ) ≤ L (j) 1 g (j) 1 (x ′ 1 γ1 , ε 1 ) if L (j) 1 < g (j) 2 (x ′ 2 γ2 , ε 2 ) < L (j) 2 L * (j) 2 if g (j) 2 (x ′ 2 γ2 , ε 2 ) ≥ L (j) 2 (2.3)
where L * (j) 1

and L * (j) 2

are qualitative measures of specific situations, L

1 and L (j) 2 are two thresholds of the selection equation. This model is a multivariate extension of the two-limit Tobit model (see e.g. [START_REF] Maddala | Estimation of Limited-dependent Variable Models under Rational Expectations[END_REF]. In addition to the two-limit selection model, it might also be useful to consider more general selection models with multiple non-response categories, such as refusals, don't know, etc., with a distinct selection equation for each category.

Real examples for potential application of the proposed model. Semiparametric MSSM has many possible real applications in economics. For example, it can be used to study the determinants of innovation behaviour or financial choices. The Community Innovation Survey collects data on the innovative characteristics of EU firms. The data include measures of innovation and related expenditures (Intramural R&D, extramural R&D, Acquisition of machinery, equipment and software and other external knowledge). MSSM could be useful to exploit this information. The selection equation could give the state "observed / non-observed" of the dependent variable y (having innovation activities) and the outcome equation would give the value of dependent variables (the amount of expenditure for each of the four innovation activities) when innovation activities are observed.

Semiparametric MSSM could also be used in clinical study when the applied researcher considers relative potency. For instance, consider a clinical study of two related drugs A and B that belong to the same class (such as two statins), with the primary goal to determine the relative potency for the two drugs. In this kind of application, it is reasonable to assume that the relative potency is determined biologically by the intrinsic nature of the two drugs, therefore the same relative potency (that is the same γj coefficients) holds for various components of the multivariate outcome measure.

Population and sample approaches

Our approach splits into two principal steps. In the first step, the idea is to use multivariate sliced inverse regression in order to get a Σ-orthogonal basis of the e.d.r. space E = Span(γ 1 , γ 2 ). In the second step, since the linear subspaces E 1 and E 2 are known (because the matrices A 1 and A 2 are assumed to be known a priori), two canonical analyses of the couples (E, E 1 ) and (E,

E 2 ) can provide bases of E ∩ E 1 = Span(γ 1 ) and E ∩ E 2 = Span(γ 2 ).

Population version

Step 1: Pooled marginal sliced inverse regression. For model (2.2), [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF] has shown that pooled marginal sliced inverse regression based on the SIR α approach, named PMS α hereafter, provides a basis denoted

B = [v 1 , v 2 ]
of the e.d.r. space E, that is Span(B) = E. The major novelty is to consider a transformation (slicing) T j (.) of y (j) with a specific slice for the missing value (MV) of y (j) . The vectors b k are the eigenvectors corresponding to the two largest eigenvalues of a Σ-symmetric matrix.

More precisely, let us give a brief overview of the PMS α method. The idea of this method is to consider the q univariate SIR α methods of each component y (j) of y on x (based on a specific slicing T j ) and to combine the corresponding M α matrices (denoted by M (j) α j ) in the following pooling:

M α,P = q j=1 w j M (j) α j , (3.1)
for positive weights w j and parameters α j ∈ [0, 1]. In the M α,P matrix, the α index stands for the vector (α 1 , . . . , α q ) and the P index stands for "pooled".

Each transformation T j categorizes each response y (j) into a new response with H j + 1 levels. We assume that the support of each y (j) is partitioned into H j fixed slices s

(j) 1 , . . . , s (j) h , . . . , s (j)
H j , plus one slice s (j) 0 for the missing value of y (j) . For j = 1 . . . , q, the matrices M (j) α j are defined as follows:

M (j) α j = (1 -α j )M (j) I Σ -1 M (j) I + α j M (j) II , with M (j) I = Var(E(x|T j (y (j) ))) = H j h=0 p (j) h (m (j) h -µ)(m (j) h -µ) ′ , M (j) II = E Var(x|T j (y (j) )) -E(Var(x|T j (y (j) ))) Σ -1 Var(x|T j (y (j) )) -E(Var(x|T j (y (j) ))) ′ = H j h=0 p (j) h V (j) h -V (j) Σ -1 V (j) h -V (j) ,
where p (j)

h = P (y (j) ∈ s (j) h ), m (j) h = E(x|y (j) ∈ s (j) h ), Var (j) h = Var(x|y (j) ∈ s (j) h )
and

V (j) = H j h=0 p (j) h V (j) h . The matrix M (j)
I is the usual matrix used in the classical SIR approach, often named SIR-I because it relies on a property of the first inverse conditional moment of x given y, while M (j)

II correspond with the SIR-II approach using information from the inverse conditional variance of x given y.

When α j = 0 (resp. α j = 1), the method used with M (j) α j is equivalent to the SIR-I (resp. SIR-II) approach for the j-th component of y.

For model (2.2), two crucial conditions for the theoretical success of SIR α and PMS α methods are the following: a linearity condition

E(v ′ x|γ ′ 1 x, γ ′ 2 x) is linear for any v, (3.2)
and a constant variance condition

Var(x|γ ′ 1 x, γ ′ 2 x) is non-random. (3.3)
Note that (3.2) is satisfied when x has an elliptically symmetric distribution and (3.3) is satisfied when x follows a multivariate normal distribution (which is an elliptical distribution). Moreover, some mild departure from the elliptical symmetry will not affect the application of SIR or MSIR, see for instance Li (1991Li ( , 1997)). Note also that low-dimensional projections from high-dimensional data are known to be able to improve the elliptical symmetry of data distribution, see for details [START_REF] Diaconis | Asymptotics of graphical projection pursuit[END_REF] or [START_REF] Hall | On almost linearity of low dimension projection from high dimensional data[END_REF]. Finally, an insightful discussion about the SIR methodology and applications can be found in [START_REF] Chen | Can SIR be as popular as multiple regression[END_REF], and most of these comments are still valid for MSIR approach.

Under conditions (3.2) and (3.3), the eigenvectors v 1 , v 2 associated with the largest two eigenvalues of Σ -1 M α,P are e.d.r. directions and span the e.d.r.

space.

Step 2a: Two canonical analysis.

Let us consider the two subspaces E k and E of R p equipped with the inner product Σ. Canonical analysis is a useful tool to find out a Σ-orthogonal basis of

E k ∩ E.
This basis is formed by the eigenvectors corresponding to the eigenvalue 1 of P E k P E , where P E k and P E are respectively the Σ-orthogonal projectors onto

E k and E.

Specifically, we have:

P E = B(B ′ ΣB) -1 B ′ Σ = BB ′ Σ and P E k = A k (A ′ k ΣA k ) -1 A ′ k Σ. It is equivalent and simpler to diagonalize P E k P E P E k which is a Σ-symmetric matrix.
Let us call b k the unique eigenvector corresponding to the eigenvalue

1 of P E k P E P E k . From Theorem 1, the eigenvector b k is colinear to γ k and is Σ-normalized: b ′ k Σb k = 1.
Step 2b: Retrieval of the direction of γk . We can derive a vector, bk , colinear to γk : bk

= A ′ k b k . This vector bk is Σ k -normalized: b′ k Σ k bk = 1.

Estimation of the directions

As mentioned in the preceding section, the directions are obtained from computations based only on covariance matrices. Substituting estimates in place of these matrices yields estimated directions.

Let {(y i , x i ), i = 1, ..., n} be a sample from the reference model (2.1). Let Σ be the empirical covariance matrix of the x i 's.

Step 1: Estimating a basis of the e.d.r. space E by PMS α method. We have to estimate the matrix M α,P . To do this using the H j + 1 slices of each component y (j) , it is straightforward to estimate the matrices M (j) I and M (j)

II by substituting empirical versions of the moments for their theoretical counterparts, and therefore to obtain the estimated matrices M (j) αj . Note that, for the choice of the slices of T j , s (j) 0 contains the cases corresponding to the missing value (MV) of y (j) . The other slices, s The choice of number H j of slices is less crucial than the choice of the smoothing parameter in nonparametric regression: in practice, we propose to choose

H j such that 2 < H j < [n * j /2],
where n * j is the number of non missing y (j) i in the sample and [a] denotes the integer part of a. For the choice of the weights w j , we can use equal weights w j = 1/q for j = 1, . . . , q if we have no a priori information on the importance of each component y (j) of y. The parameters α j are individually chosen for each matrix M (j)

α j , and we propose to use the method based on the test approach of [START_REF] Saracco | Pooled Slicing methods versus Slicing methods[END_REF], which does note require the estimation of the link functions. Therefore we obtain the estimated matrix:

Mα,P = 1 q q j=1 M (j) αj .
(3.4)

The two estimated e.d.r. directions, v1 and v2 , are then the eigenvectors corresponding to the two largest eigenvalues of Σ-1 Mα,P . These vectors form a Σ-orthonormal system. Let Ê = Span( B) where

B = [v 1 , v2 ].
Step 2a: Estimating the direction of γ k , k = 1, 2. We obtain these directions by canonical analyses of ( Ê, E 1 ) and ( Ê, E 2 ). Thefore, the estimate of the direction of γ k is the eigenvector bk corresponding to the major eigenvalue of the Σ-symmetric matrix PE k P Ê PE k , where

P Ê = B( B′ Σ B) -1 B′ Σ = B B′ Σ and PE k = A k (A ′ k ΣA k ) -1 A k Σ.
Step 2b: Estimating the direction of γk , k = 1, 2. The estimates of the direction of γk are then given by bk = A ′ k bk . Remarks: -In order to obtain an estimate of the entire vector γk (and not only of its direction), we can normalize this vector in the Σ k metric, and impose the sign of a non null component of γk .

-For the two-limit model (2.3), there must be one slice for each kind of missing y value. The other slices are built, splitting the other cases in the usual way.

We will study the asymptotic properties of the estimators b1 and b2 in the next section. First, however, we will discuss a topic of practical concern connected with the estimation process: the estimation of the link functions of the model (2.1).

Rough approximation of the link functions g (j)

1 and estimation of the state of y probabilities. Let us simplify the reference model by assuming an additive error component: for j = 1, . . . , q, y (j) = g

(j) 1 (x ′ 1 γ1 ) + ε (j) 1 if g (j) 2 (x ′ 2 γ2 ) + ε (j) 2 > 0 0 otherwise. with E(ε (j) 1 ) = E(ε (j)
2 ) = 0. A rough approximation of the jth-observation link function, g

(j)
1 , may be nonparametrically obtained by kernel or spline methods. [START_REF] Eubank | Spline smoothing and nonparametric regression[END_REF] and [START_REF] Haerdle | Applied nonparametric regression[END_REF] give an operational description of these tools.

We may, for instance, build a naive Nadaraya-Watson kernel estimate from the subsample of cases where y (j) is non missing by regressing y (j) on x′ 1 b1 . This estimator is generally a biased estimator of g (j)

1 since E(ε

(j) 1 |x ′ 1 γ1 , g 2 (x ′ 2 γ2 )+ε (j) 2 > 0) is non null.
Let us now examine the case of the selection link functions, g (j)

2 . What is interesting is to estimate the probability of the state of y (j) . In order to describe the state of y (j) , let us introduce the qualitative variable t (j) for the one-limit selection model (2.1) (resp. for the two-limits selection model (2.3)):

t (j) = 1 if y (j) is observed 0 otherwise   resp. t (j) =      0 if y (j) = L * 1 1 if y (j) is observed 2 if y (j) = L * 2    .
From each sample {(t

(j)
i , r i ), i = 1, . . . , n} where r i = x′ 2i b2 , we can obtain a naïve Nadaraya-Watson estimate of the probability P (t (j) = t|r = x′ 2 b2 ) by:

p(j) n (t|r) = n i=1 K r-r i νn n l=1 K r-r l νn I[t (j) i = t], (3.5)
where K is a kernel function and ν n is the bandwidth which may be chosen by cross validation.

Asymptotic theory

In the sequel, the notation X n -→ d X means that X n converges in distribution to X as n → ∞. Let D 1 ⊗ D 2 denote the Kronecker product of the matrices D 1 and D 2 (see [START_REF] Tyler | Asymptotic inference for eigenvectors[END_REF] for some useful properties of the Kronecker product). From now on, for each s × s matrix s1) , d (21) , d (22) , . . . , d (ss) ′ be the s 2 -dimensional column vector of all elements of D.

D = d (jk) , let vec(D) = d (11) , . . . , d (
The necessary assumptions are gathered together below for easy reference.

(A1) {(y i , x i ), i = 1, . . . , n} is a sample of independent observations from model (2.1).

(A2) The supports of each component y (j) (when observed) of y are partitioned into H j fixed slices s

(j) 1 , . . . , s (j)
h , . . . , s

H j such that p (j) h = 0, with a special slice s (j) 0 for the missing y (j) .

(A3) The covariance matrix Σ is positive definite.

(A4) The two largest eigenvalues of Σ -1 M α,P satisfy λ 1 ≥ λ 2 > λ 3 ≥ 0. Proof. Classical asymptotic theory gives us: Σ = Σ + O p (n -1/2 ). By the asymptotic theory of PMS α (see [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF], we get B = B + O p (n -1/2 ). Thus,

P Ê = P E + O p (n -1/2 ). (4.1)
From the identifiability conditions, rank

(A ′ k ΣA k ) = p k . Since A ′ k ΣA k = A ′ k ΣA k + O p (n -1/2 ), we get (A ′ k ΣA k ) -1 = (A ′ k ΣA k ) -1 + O p (n -1/2 ) and PE k = P E k + O p (n -1/2 ), j = 1, 2. (4.2) Combining (4.1) with (4.2) yields PE k P Ê PE k = P E k P E P E k +O p (n -1/2 ), k = 1, 2.
Consequently, the eigenvector of PE k P Ê PE k corresponding to the major eigenvalue converges at the same rate to the corresponding eigenvector for A3) and (A4), we have:

P E k P E P E k : bk = b k + O p (n -1/2 ), k = 1, 2.
for k = 1, 2, √ n( bk -bk ) -→ d N (0, A ′ k G k C * G ′ k A k ),
where the expression of G k is given in (4.3) and the expression of C * can be found in [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF].

Proof. The proof is divided into three steps.

Step 1: Asymptotic distribution of the Canonical analysis matrix.

Let us consider the decomposition √ n(

PE k P Ê PE k -P E k P E P E k ) = √ n( PE k -P E k ) P Ê ( PE k -P E k ) + √ n( PE k -P E k ) P Ê P E k + √ nP E k P Ê ( PE k -P E k ) + √ n(P E k P Ê P E k -P E k P E P E k ).
The first term of the right hand side is O p (n -1/2 ). Thus,

√ n vec( PE k P Ê PE k ) -vec(P E k P E P E k )
has the same asymptotic distribution as the last three terms of the decomposition. These terms can be written as follows:

([P ′ E k P ′ Ê ⊗I p ]+[I p ⊗P E k P Ê ]) √ n[vec( PE k )-vec(P E k )]+[P ′ E k ⊗P E k ] √ n[vec( P Ê )-vec(P E )].
We prove in the Appendix that √ n vec( PE k )vec(P E k ) has the same asymp-

totic distribution as N k √ n[vec( Σ) -vec(Σ)]
where N k is defined in (6.1). Moreover, from [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF], since Σ-1 Mα,P converges in probability to Σ -1 M α,P , we have, with a probability converging to 1, for n sufficiently large: || Σ-1 Mα,P -

Σ -1 M α,P || ≤ λ 2 /2
, where λ 2 is the second major eigenvalue of Σ -1 M α,P . Then we can apply the Lemma 4.1 of [START_REF] Tyler | Asymptotic inference for eigenvectors[END_REF], and we obtain the asymptotic distribution of the eigenprojector on the estimated e.d.r. space:

√ n vec( P Ê ) -vec(P E )
has the same asymptotic distribution as

C w √ n vec( Σ-1 Mα,P ) -vec(Σ -1 M α,P )
where

C w = -λ∈w [(M α,P Σ -1 -λI p ) + ⊗ P λ + P ′ λ ⊗ (Σ -1 M α,P -λI p ) + ], with w = {λ 1 , λ 2 }.
Finally, the asymptotic distribution of √ n(

PE k P Ê PE k -P E k P E P E k ) is then the same as Â0 √ n vec( Σ-1 Mα,P ) vec( Σ) - vec(Σ -1 M α,P ) vec(Σ) , where Â0 = A 0 1 | Â0 2 , with A 0 1 = (P ′ E k ⊗P E k )C w and Â0 2 = [P E k P ′ Ê ⊗ I p ] + [I p ⊗ P E k P Ê ] N k . Moreover, it is easy to show that Â0 -→ P A 0 where A 0 = A 0 1 | A 0 2 with A 0 2 = [P E k P ′ E ⊗ I p ] + [I p ⊗ P E k P E ] N k .
Step 2: Asymptotic distribution of the major eigenvector.

Remembering that bk (resp. b k ) is the eigenvector corresponding to the major eigenvalue of PE k P Ê PE k (resp. P E k P E P E k ), we apply Lemma 2 of [START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF].

First, we need to specify the asymptotic distribution of

√ n vec( PE k P Ê PE k ) vec( Σ) - vec(P E k P E P E k ) vec(Σ) .
From step 1, this vector has the same asymptotic distribution as

B0 √ n vec( Σ-1 Mα,P ) vec( Σ) - vec(Σ -1 M α,P ) vec(Σ) , where B0 = A 0 1 Â0 2 0 p 2 ,p 2 I p 2 . Since Â0 2 -→ P A 0 2
, we get B0 -→ P B 0 where

B 0 = A 0 1 A 0 2 0 p 2 ,p 2 I p 2 .
Moreover, from an application of the Delta method, [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF] shows that

√ n vec( Σ-1 Mα,P ) vec( Σ) - vec(Σ -1 M α,P ) vec(Σ) -→ d Φ * = vec(Φ) vec(Φ Σ ) ∼ N (0, C * ),
The expression of C * can be found in [START_REF] Saracco | Asymptotics for pooled marginal slicing estimator based on SIR α approach[END_REF].

Thus we obtain:

√ n vec( PE k P Ê PE k ) vec( Σ) - vec(P E k P E P E k ) vec(Σ) -→ d B 0 Φ * ,
where B 0 Φ * ∼ N (0, B 0 C * B 0′ ). We can now apply Lemma 2 of [START_REF] Saracco | An asymptotic theory for Sliced Inverse Regression[END_REF],

and we get:

√ n( bk -b k ) -→ d R k , where R k = [b ′ k ⊗ (P E k P E P E k -I p ) + ]B 0 vec(Φ) vec(Φ Σ ) -1 2 (b ′ k Φ Σ b k )b k . Tedious but
simple computations give us for R k a multivariate normal distribution with mean zero and covariance matrix G k C * G ′ k where the matrix G k is:

{b ′ k ⊗ (P E k P E P E k -I p ) + }(P ′ E k ⊗ P E k )C w | {b ′ k ⊗ (P E k P E P E k -I p ) + }((P E k P E ) ′ ⊗ I p + I p ⊗ P E k P E )N k -1 2 b k (b ′ k ⊗ b ′ k ) .
(4.3)

Step 3: Asymptotic distribution of bk .

Finally, since bk = A ′ k bk and bk = A ′ k b k , we get:

√ n( bk -bk ) -→ d Rj = A ′ k R k where Rk ∼ N (0, A ′ k G k C * G ′ k A k ). Remark.
From a theoretical point of view, the asymptotic covariances of these two estimators can be estimated by replacing the theoretical terms by their empirical √ n-consistent counterparts. The corresponding estimated asymptotics matrices converge to the true ones at rate √ n. From a computational point of view, it is tedious to obtain these estimators of the asymptotic covariances. Nevertheless, we can easily compute bootstrap estimators of these matrices which are very close to the true matrices (obtained by Monte-Carlo method). We illustrate this point in Section 5.1 on a simulated example.

Simulation results

In order to evaluate the numerical performance of the proposed method, a simulation study was carried out. Following [START_REF] Duan | Slicing regression: a link-free regression method[END_REF], we measure the quality of the estimate bk of the direction of γk by:

cos 2 bk , γk = b′ k Σ k γk 2 ( b′ k Σ k bk )(γ ′ k Σ k γk ) where Σ k = A ′ k ΣA k .
The closer the squared cosine is to one, the better the estimation.

We generate simulated data from the semiparametric multivariate (q = 2) model (2.1) with

           g (1) 1 (x ′ 1 γ1 , ε (1) 1 ) = exp(x ′ 1 γ1 ) + ε (1) 1 g (1) 2 (x ′ 2 γ2 , ε (1)
2

) = x′ 2 γ2 + ε (1) 2 g (2) 1 (x ′ 1 γ1 , ε (2) 1 ) = (x ′ 1 γ1 ) 3 + 3(x ′ 1 γ1 ) + ε (2) 1 g (2) 2 (x ′ 2 γ2 , ε (2)
2 ) = (x ′ 2 γ2 ) 2 + ε

(2) 2

(5.1)

where x follows a p-dimensional standardized normal distribution, and x1 (resp.

x2 ) is the (p -1)-dimensional vector corresponding to the first (resp. last) (p -1) coordinates of x. The error term ε = (ǫ

(1) 1 , ǫ (1)
2 , ǫ

(2)

1 , ǫ (2)
2 ) ′ is normally distributed: ǫ ∼ N 4 (µ ǫ , Σ ǫ ). Two designs of the covariance of ǫ will be considered:

Σ I ǫ =       1 ρ 0 0 ρ 1 0 0 0 0 1 ρ 0 0 ρ 1       and Σ II ǫ =       1 ρ ρ ρ ρ 1 ρ ρ ρ ρ 1 ρ ρ ρ ρ 1      
, with different values of ρ (0.1, 0.5 and 0.9). In the matrix Σ I ǫ the error term associated with the two components y (1) and y (2) are assumed to be independent, which is not the case with the covariance matrix Σ II ǫ . Note that we never consider the most favourable case with an independent error term between the observed equation and the selection equation. To control the number of non observed values for the y (j) 's component, we use two different values of µ ǫ : in order to obtain around 25% (resp. 50%) of non observed values for y (1) and y (2) , we choose µ ǫ = (0, 1.5, 0, -0.5) (resp. µ ǫ = (0, 0, 0, -2)). For the slope parameters, we take γ1 = (1, 1, -1, -1, 0, . . . , 0) ′ and γ2 = (0, . . . , 0, 1, -1, 1, -1) ′ .

To study the performance of the proposed method, we consider different sample sizes (n = 100, 200 and 300), various dimensions of the explanatory variable (p = 5, 10), the two different choices of covariance matrix (Σ I ǫ and Σ II ǫ ), and two levels L of non observed values for y (j) (25% and 50%). The number of slices in the PMS α method, H j , is specified to be H j = max n * j , p where n * j is the number of observed y

i 's in the sample. In the next subsection, we apply our approach to a simulated sample. Then, in the last subsection, we comment on the complete simulation study. Simulations were performed with R. All of the source codes are available from the authors by e-mail.

Simulated example

In this susbsection, we consider a simulated sample of n = 100 data points from the previous model for p = 5, Σ ǫ = Σ II ǫ , ρ = 0.5 and L = 25%. On the left hand side of Figure 5.1, we can observe the plots of the response variables y (1) and y (2) versus the true "observation" index x′ 1 γ1 . Let us introduce the two variables y

(1) * = g (1) 2 (x ′ 2 γ2 , ε (1)
2 ) and y

(2) * = g (2) 2 (x ′ 2 γ2 , ε (2)
2 ), which are called in the literature latent variables (since in practice the values of these variables are never available in the sample). On the right hand side of Figure 5.1, we plot these latent variables y (j) * versus the true "selection" index x′ 2 γ2 . The horizontal line allows us to determine for which observations the y (j) i 's values will be non observed in the left hand side graphics. The directions of γ1 and γ2 are then estimated and we get b1 = (-0.483, -0.565, 0.447, 0.497) ′ and b2 = (-0.613, 0.539, -0.350, 0.459) ′ . The corresponding squared cosines are respectively equal to 0.993 and 0.962. Note that b1 (resp. b2 ) gives nearly the same direction as γ1 (resp. γ2 ). Moreover, we compute the quality of the estimation Ê of the e.d.r. space E using Trace(P E P Ê )/2 which is equal to 0.886 for this simulated sample. Even if this subspace is relatively poorly estimated compared with the quality of each estimated direction, the second step (which takes into account additional information) ensures that we recover the good directions of the observation and selection slope vectors. We estimate the asymptotic covariance matrices, denoted by V ( b1 ) and V ( b2 ), with the bootstrap method (with 500 replications):

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !4 !2 0
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !4 !2 0 2 4 !2 0 2 4 6 true selection index y * ( (1) ) ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !4 !2 0 2 4 !40 0
! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !4 !2 0
V ( b1) = 10 -3 0 B B B @
4.17 -1.93 1.74 -0.90 4.00 0.24 4.16 2.84 -0.36 17.1 These matrices are very close to the "true" asymptotic covariance matrices,V ( b1 ) and V ( b2 ) (not given here), calculated via the Monte Carlo approach. Note that the variance terms in V ( b2 ) are greater than those obtained in V ( b1 ), because of the low level L (=25%) of non observed values for y (j) .

1 C C C A , V ( b2) = 10 -2 0 B B B @ 4.
In Figure 5.2, we represent on the left hand side the plots of the response variable y (j) versus the estimated "observation" index x′ 1 b1 . Note that, since we have b1 ≃ -γ 1 /||γ 1 || (resp. b2 ≃ -γ 2 /||γ 2 ||), the scatterplots of Figures 5.1 and 5.2 (left hand side) do not have the same orientation. We add on these plots the Nadaraya-Watson estimate of the observation link functions. On the right hand side, we plot the t (j) 's values versus the estimated "selection" index x′ 2 b2 , and we also plot the Nadaraya-Watson estimate of the probability to observe y (j) , based on the equation (3.5).

! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !3 !2 !1 0 1 2 3 0 10 30 50 estimated observation index y ( (1) ) !3 !2 !1 0 1 2 3 0 10 30 50 estimated observation index y ( (1) ) ! observed values non observed values ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! !2 !1 0 1 2 0.0 0.4 0.8 estimated selection index t ( (1) ) ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !2 !1 0 1 2 !40 0 20 40 estimated observation index y ( (2) ) !2 !1 0 1 2 !40 0 20 40 estimated observation index y ( (2) ) ! observed values non observed values ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !! ! ! ! ! ! ! ! ! ! ! ! !2 !1 0 1 2 0.0 0.4 0.8 estimated selection index t ( (2) )
Figure 5.2: Kernel estimate of the observation link functions (left hand side) and Nadaraya-Watson estimate of the probability of t (j) = 1 (that is y (j) observed)

Results of the simulation study

In our study, we consider combinations of the following simulation parameters: the level L of non observed values for y (j) (25% or 50%), the form of the error covariance matrix Σ ǫ (Σ I ǫ or Σ II ǫ with ρ = 0.1, 0.5 or 0.9, and the dimension p of the covariable (p = 5 or 10). We also take into account various sample sizes n = 100, 200 or 300.

For each combination, N = 500 samples were generated. For each sample l = 1, . . . , N , the directions of the slope vectors γ1 and γ2 were estimated and we get bl 1 and bl 2 . Then, we evaluate the corresponding values of the quality measure: c l k = cos 2 ( bl k , γk ) for k = 1, 2 and l = 1, . . . , N . We show the results via the boxplots of these squared cosines for different combinations. When p = 5 and Σ ǫ = Σ I ǫ (resp. Σ ǫ = Σ II ǫ ), Figure 5.3 (resp. Figure 5.4) gives the boxplots for γ1 and γ2 denoted by G1 and G2 in the graphics, for the different values of ρ, L and n. Figure 5.5 shows the boxplots when p = 10, Σ ǫ = Σ II ǫ and n = 300, for various ρ; note that the vertical scale in this figure goes from 0.75 to 1 (contrary to the previous one which goes from 0.4 to 1).

From Figures 5.3, 5.4 and 5.5, we can see that the results with these simulated data are very good. More precisely, one can observe that:

-The estimations of the γ1 and γ2 's directions are good since almost all boxplots of the squared cosines are concentrated in the interval [0.9, 1].

-The form of the covariance matrix of the error term ε and the value of the parameter ρ do not seem to have any influence on the quality of the estimates.

-The level L of the non observed values for the y (j) 's has only a slight influence on the quality of the estimation of the selection slope vectors γ2 , especially in terms of spread of the squared cosine values. When this level is low (L = 25%), there is less information on the selection part of the model so the quality of the γ2 estimates is slightly lower than when this level is larger (L = 50%). On the other hand, not surprisingly, there is an opposite behavior for the estimates of the observation slope parameter γ1 since there is less information on the observation part of the model when L is large.

-The sample size n has a quite predictable influence of the quality of the estimates: the largest is the sample size, the greatest are the squared cosines. When n = 200 or 300, the quality of the two estimated directions are very good.

-Dimension p of the explanatory variable x does not seem to have any effect on the quality of the estimates.

Simulation with a non-normal distributed covariable x

In order to investigate the robustness of the method when x does not fol-

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.1, 25% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.5, 25% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.9, 25% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.1, 50% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.5, 50% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0 ρ ρ = = 0.9, 50% for µ ǫ to obtain around 25% (resp. 50%) of non observed values for the y (1) 's and the y (2) 's. Moreover we take γ1 = (1, 1, -1, -1) ′ /2 and γ2 = (1, -1, 1, -1) ′ /2 for the observation and selection slope parameters.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ρ ρ = = 0.1, 25% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.5, 25% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.9, 25% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.1, 50% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0

ρ ρ = = 0.5, 50% G 1 100 G 2 100 G 1 200 G 2 200 G 1 300 G 2 300 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
0.4 0.5 0.6 0.7 0.8 0.9 1.0 ρ ρ = = 0.9, 50% ! ! = = 0.9 G1, 25% G2, 25% G1, 50% G2, 50% (not detailed here) were observed for the two other distributions.

Comparison with a parametric approach

We compare in this simulation the parametric Tobit II model (implemented in [START_REF] Henningsen | sampleSelection -A package for the free Statistical Software[END_REF]) with our semiparametric approach using various error distributions and various selection and observations link function shapes. We consider here two models, (M 1) and (M 2), from the sample models defined in (2.1) with q = 1:

(M 1) :

g 1 (x ′ 1 γ1 , ǫ 1 ) = x′ 1 γ1 + ǫ 1 g 2 (x ′ 2 γ2 , ǫ 2 ) = x′ 2 γ2 + ǫ 2
and (M 2) :

g 1 (x ′ 1 γ1 , ǫ 1 ) = exp (x ′ 1 γ1 ) + ǫ 1 g 2 (x ′ 2 γ2 , ǫ 2 ) = exp (x ′ 2 γ2 ) + ǫ 2
Model (M 1) is in favour of the parametric approach with linear link functions, whereas model (M 2) has non-linear link functions. For these two models, the error term ǫ = (ǫ 1 , ǫ 2 ) is normally distributed as in the previous simulation study,

x follows a five-dimensional standardized normal distribution, x1 (resp. x2 ) is the 4-dimensional vector corresponding to the first (resp. last) four coordinates of x.

To control the level L of non observed values for y, we used different values of µ ǫ .

For the slope parameters , we took γ1 = (1, 1, -1, -1)

′ and γ2 = (1, -1, 1, -1) ′ .

We present in Figure 5.7 only the results for n = 200, ρ = 0.9 and L = 50%, over N = 500 replicated samples. As is to be expected, the Tobit II approach three referees for their valuable comments and constructive suggestions. They thank Jean Belin, economics researcher of the GREThA laboratory of Bordeaux 4

University for numerous discussions and future work on applications in economics area of the proposed semiparametric multivariate sample selection model. K p is the p 2 × p 2 commutation matrix (see [START_REF] Magnus | The commutation matrix: some properties and applications[END_REF].

Appendix: Asymptotic distribution of PE k Let PE k = A k (A ′ k ΣA k ) -1 A ′ k Σ ( resp. P E k = A k (A ′ k ΣA k ) -1 A ′ k Σ)
We obtain the asymptotic distribution of PE k through the following three steps.

Step 1. Let f 1 : R p 2 -→ R p k p+p 2 k be defined by

f 1 (vec(M )) = (I p ⊗ A ′ k )vec(M ) (A ′ k ⊗ A ′ k )vec(M )
.

Then, from the Delta method, we get:

√ n vec(A ′ k Σ) vec(A ′ k ΣA k ) - vec(A ′ k Σ) vec(A ′ k ΣA k ) -→ d U 1k ,
where U 1k ∼ N (0, C 1k ) with

C 1k = I p ⊗ A ′ k A ′ k ⊗ A ′ k C Σ I p ⊗ A k A k ⊗ A k .
Step 2. From the following first order approximation:

√ n((A ′ k ΣA k ) -1 -(A ′ k ΣA k ) -1 ) . = -(A ′ k ΣA k ) -1 √ n(A ′ k ΣA k -A ′ k ΣA k ) (A ′ k ΣA k ) -1 ,
we derive:

√ n vec(A ′ k Σ) vec((A ′ k ΣA k ) -1 ) - vec(A ′ k Σ) vec((A ′ k ΣA k ) -1 ) -→ d U 2k = S k U 1k
where U 2k ∼ N (0, C 2k ) with C 2k = S k C 1k S ′ k and S k =

I p k p 0 p k p+p 2 k 0 p 2 k +p k p -(A ′ k ΣA k ) -1 ⊗ (A ′ k ΣA k ) -1 .
Step 3. Let us introduce the function f 2 : R p k p+p 2 k -→ R p 2 defined by f 2 vec(M 1 ) vec(M 2 ) = vec(A k M 2 M 1 ). Then from a second application of the Delta method, we derive:

vec( √ n[ Pj -P E k ]) -→ d U k , where U k ∼ N (0, C U k ) with C U k = N k C Σ N ′ k and N k = I p ⊗ A k (A ′ k ΣA k ) -1 A ′ k -P ′ E k ⊗ A k (A ′ k ΣA k ) -1 A ′ k = (I p -P ′ E k ) ⊗ [A k (A ′ k ΣA k ) -1 A ′ k ]
(6.1)

  = 1, ..., H j , are made by splitting the range of the non-missing values of the jth component of y into slices of nearly equal weight.

4. 1 .

 1 Convergence in probability of the estimated directions Theorem 2. Under conditions given in (3.2) and (3.3), and under assumptions (A1), (A2) and (A3), we have bk = bk + O p (n -1/2 ), with the vector bk colinear to γk , for k = 1, 2,.

  Finally, since bk = A ′ k bk and bk = A ′ k b k , we conclude that: bk = bk + O p (n -1/2 ), k = 1, 2. From Theorem 1, we have bk colinear to γk . 4.2. Asymptotic distribution of bk , k = 1, 2 Theorem 3. Under conditions (3.2) and (3.3), and under assumptions (A1), (A2), (

Figure 5 . 1 :

 51 Figure 5.1: Plots of y (j) versus the true "observation" index x′ 1 γ1 (on the left) and plots of the latent variables y (j) * versus the true "selection" index x′ 2 γ2 (on the right).

Figure 5 . 3 :

 53 Figure 5.3: Boxplots of the squared cosines when Σ ǫ = Σ I ǫ and p = 5

Figure 5 . 4 :

 54 Figure 5.4: Boxplots of the squared cosines when Σ ǫ = Σ 2 ǫ and p = 5

Figure 5 .

 5 Figure5.6 reports the results of this simulation study via the boxplots of these squared cosines for the discrete rectangular distribution. One can see that the estimations of the directions of the slopes for the selection equations and the outcome equations are quite good, even for a discrete x that does not follow an elliptically symmetric distribution. Note that, as in the multivariate normal case (in the previous subsection), we can observe the same influence of the level L on

Figure 5 . 5 :

 55 Figure 5.5: Boxplots of the squared cosines when Σ ǫ = Σ 2 ǫ , n = 300 and p = 10

  be the Σ (resp. Σ) orthogonal projector onto the linear subspace E k spanned by the columns of A k .For an elliptically distributed x, with covariance matrix Σ and kurtosis parameter κ,[START_REF] Tyler | Asymptotic inference for eigenvectors[END_REF] gave the following asymptotic distribution:√ n( Σ -Σ) -→ d Φ Σwhere vec(Φ Σ ) ∼ N (0, C Σ ) and C Σ = (1+κ)(I p 2 +K p )(Σ⊗Σ)+κvec(Σ)[vec(Σ)] ′ .

Acknowledgment

The authors are very grateful to the editor, the associate editor and the

We have also considered various combinations of the simulation parameters: the level L of non observed values for y (j) (25% or 50%), the error term correlation ρ = 0.1, 0.5 or 0.9. In any case, we observed very similar results. Moreover, we compared the two approaches with the linear model (M 1) when the error term is non normal distributed. The Tobit II method appears to be robust to mild violations of the normality assumption like our approach (which does not rely on this kind of assumption).

Concluding remarks

In this paper, we proposed a new semi-parametric estimation method for a multivariate sample selection model (MSSM). As pointed out previously, our semi-parametric approach has the main advantage of being link-and distributionfree. The proposed geometric approach to the estimation of the slope vectors in

0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.4 0.5 0.6 0.7 0.8 0.9 1.0 ρ ρ = = 0.9, 50% the outcome equation and in the selection models has also the advantage of dealing symmetrically with both slope vectors. From a theoretical point of view, the convergence in probability at root n rate and the asymptotic normality of the slope estimators have been proved. This estimation method is numerically very fast since it is based on only a few matrix calculus and eigen-decompositions and does not demand any time-consuming iterative computations. Note that the corresponding algorithm is easy to implement. The R source code is available from the authors. Another interesting aspect is that this method does not require starting values. From a practical point of view, the simulation study has highlighted a good behaviour of the estimation method even for non-elliptical distribution of the covariate. Moreover a real economic application is currently under investigation. Finally a direction which would probably be interesting to investigate would be to develop another two-step semi-parametric estimation methods. In a first step, we could take into account the MSIR estimator of the selection slope parameter since the selection probability only depends on the index x2 γ2 . Then in a second step, we will incorporate the additional information in order to get the observation slope vector from the entire e.d.r. space.


Sabine Beuermann

Simon Harrisson

Robin A Hutchinson

Tanja Junkers

Gregory T Russell

Update and critical reanalysis of IUPAC benchmark propagation rate coefficient data

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Between 1995 and 2017, the IUPAC Polymerization Kinetics Subcommittee published a series of papers [1][2][3][4][5][6][7][8] providing benchmark values for the activation energies (EA) and preexponential factors (A) of propagation rate coefficients, kp, of 13 common monomers in radical polymerization. These rate coefficients were measured using the pulsed laser polymerization (PLP) technique, 9 in which a laser is used to generate periodic bursts of radicals in a solution of monomer maintained at a controlled temperature for a sufficient time to convert a small fraction of the monomer (generally <2%) to polymer. A substantial fraction of the resulting polymer chains is initiated by one burst and terminated by the following burst, which are separated by a time, t0, leading to a characteristically shaped polymer molar mass distribution (MMD) which depends on the propagation rate coefficient of the monomer under investigation. Specifically, a distinctive peak in the MMD corresponding to the chain length DP0 is directly proportional to kp, t0, and monomer concentration [M] according to the relation DP0 = kp[M]t0. With t0 and [M] precisely known, the accuracy of the technique is directly dependent on the measurement of DP0, typically determined from analysis of the polymer MMD measured by size exclusion chromatography (SEC). Even with careful SEC analysis utilizing the principal of universal calibration verified by multi-detector analysis, the principal source of uncertainty in the estimation of kp arises from the polymer analysis, as detailed in the previous studies. [1][2][3][4][5][6][7][8] The benchmark values were obtained by pooling the results of multiple laboratories (from 2 to 9 depending on the study) and verifying that certain experimental conditions (e.g. invariability with respect to pulse rate, presence of at least one overtone, etc.) were fulfilled. The resulting dataset has been highly useful to the polymerization community and underpin many further kinetic studies and simulations, which is demonstrated by their remarkable number of citations. These and other selected kp data have recently been collected in a machine-readable database. 10 Since the publication of the benchmark dataset for styrene more than 25 years ago, many further PLP studies have been carried out according to the IUPAC guidelines. These provide additional, independent data that can be used to refine the benchmark parameter estimates. Additionally, by pooling the results of 83 individual PLP studies on the 13 monomers, in this contribution we are able to provide an estimate of the interlaboratory variation resulting from systematic errors that are constant within a single study. Taking this variation into account, we arrive at larger but more realistic estimates of uncertainty than those originally published. It is important to note that while the original parameter estimates differ slightly from the revised values presented here, they remain within the margin of uncertainty for the updated estimates.

Results and Discussion

IUPAC-benchmarked activation energies and pre-exponential factors are available for the monomers styrene (STY), 1 methyl methacrylate (MMA), 2 ethyl methacrylate (EMA), 3 n-butyl methacrylate (BMA), 3 n-dodecyl methacrylate (DMA), 3 cyclohexyl methacrylate (CHMA), 4 glycidyl methacrylate (GMA), 4 benzyl methacrylate (BnMA), 4 isobornyl methacrylate (iBoMA), 4 methacrylic acid (MAA), 5 n-butyl acrylate (BA), 6 methyl acrylate (MA), 7 and vinyl acetate (VAc). 8 These values are for the bulk monomer in all cases except for that of MAA, which is given for a 15% aqueous solution in water at natural pH. The parameters, as well as the calculated kp at 25°C, are collated in Table 1, together with the estimated standard error in each parameter reported in the original studies.

Table 1. Comparison of originally published IUPAC benchmark values and re-analysed values (using expanded dataset) for Arrhenius parameters of propagation rate coefficients of vinyl monomers. Numbers in parentheses represent the standard error in the final digits: e.g. 32.5(3) represents 32.5 ± 0.3. Significant (>1σ) differences between original and revised values are highlighted in bold.

Monomer

Original A notable aspect of the IUPAC benchmarking studies was the care that was taken to provide estimates not only of EA and A, but also of their uncertainties. These were presented as 95% joint confidence regions (JCRs): an identically constructed study would be expected to produce parameter estimates within these regions 95% of the time, assuming that the experimental errors are independent and identically distributed. In all cases, the JCRs were highly correlated -the error in A depended greatly on the error in EA, leading to elongated, banana-shaped JCRs.

IUPAC Benchmark Values a Revised IUPAC Benchmark Values b A (L•mol -1 •s -1 ) EA (kJ•mol -1 ) kp at 25°C (L•mol -1 •s -1 ) T (°C) N c A (L•mol -1 •s -1 ) EA (kJ•mol -1 ) kp at 25°C (L•mol -1 •s -1 ) T (°C) N c STY 1 10 7.
The correlation between the errors in A and EA is due to nature of the Arrhenius relationship (eq. 1)

𝑘 p = 𝐴 𝑒 -𝐸 A R𝑇 (1)
Experimentally, A is determined by extrapolating experimental kp vs T data to infinite temperature. Thus, a small variation in EA will result in a large variation in the extrapolated value of A. The resulting JCRs can be difficult to compare, as the uncertainty in A is dominated by the uncertainty in EA. This correlation can be reduced by modifying the Arrhenius relationship as follows (eq. 2):

𝑘 p = 𝑘 p0 𝑒 - 𝐸 A R ( 1 𝑇 - 1 𝑇 0 ) (2)
In eq. 2, kp0 is the kp at a reference temperature T0, chosen to be within the range of experimentally accessible temperatures. This corresponds to a simple change of variables, and the preexponential factor A can be obtained by setting 1/T = 0. By appropriate choice of T0, the correlation between kp0 and EA can be greatly reduced, or even eliminated. As a result, uncertainties in the parameters of the Arrhenius relationship can be presented concisely as follows (eq. 3):

𝑘 p = (𝑘 p0 ± 𝜎 𝑘 p0 ) 𝑒 - (𝐸 A ±𝜎 𝐸 A ) R ( 1 𝑇 - 1 𝑇 0 ) (3)
where σkp0 and σEA represent the uncertainties in kp0 and EA, respectively, and T0 is the temperature at which these uncertainties are uncorrelated. ‡ The uncertainty in the preexponential factor is then obtained from the propagation of errors as (eq. 4)

𝜎 𝐴 𝐴 = √( 𝜎 𝑘 p0 𝑘 p0 ) 2 + ( 𝜎 𝐸 A R𝑇 0 ) 2 (4)
In the remainder of this paper, this representation of the Arrhenius relationship is used. The reference temperature is set at 25°C (298.15 K), a temperature which falls within the experimental datasets of all monomers under consideration. In this way, the activation energies and rate coefficients at 25°C can readily be compared.

Data treatment

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 3

Please do not adjust margins Please do not adjust margins

In establishing a set of benchmark values, the IUPAC group first needed to establish a reliable set of data points. This was done by including only data that met a set of experimental conditions including duplicate experimental runs, the presence of at least one overtone peak in the molar mass distribution, and invariance of the results with respect to changes in the radical concentration, pulse repetition rate and duration of irradiation.

The temperature range was limited to less than 90°C in the case of methacrylates [2][3][4] to avoid interference from depropagation reactions, and less than 20°C 6 or 60°C 7 in the case of acrylates to avoid interference from backbiting. 5,7,11 These curated datasets were published in the original articles, [1][2][3][4][5][6][7][8] and are reproduced in the Supporting Information.

Once the dataset had been established, further issues were encountered in analysing the data:

• Only a limited number of laboratories were suitably equipped to carry out the experiment, so relatively few laboratories participated in each study.

• The number of data points provided by each laboratory was not constant, so that there was a risk that a single laboratory that provided a large number of data points would dominate the dataset.

• Different laboratories carried out experiments over different temperature ranges.

In order to mitigate these problems, when one laboratory dominated the dataset, results from other laboratories were given extra weight. 1,3,4 Likewise, when results at high or low temperature were dominated by a single laboratory, the temperature range was restricted to exclude those results. 3,4 The remaining weighted results were pooled, and either fitted directly to eq. 1 (using nonlinear least squares fitting), or transformed by taking the natural logarithm of kp and fitting a straight line as a function of 1/T. As the kp is determined by analysis of a MMD, usually obtained by SEC in which the elution volume of a polymer is approximately proportional to the logarithm of its molecular weight, errors in ln kp are constant, and this transformation of the data does not bias the parameter estimates. 2,12 While these strategies were quite effective at reducing the impact of a single laboratory on the combined datasets, they have some disadvantages. The weighting procedure is somewhat arbitrary, while restricting the temperature range can involve discarding a significant quantity of data. The most important issue, however, is the assumption, implicit in the procedure of fitting a single line to the combined dataset, that the experimental errors in each point are uncorrelated. In practice, this is not the case.

To take styrene as an example, the IUPAC paper 1 contained data from 8 studies, of which 4 provided values at a single temperature (25°C), [13][14][15][16] and 4 provided multiple kp values across a range of temperatures. 12,[17][18][19] As shown in Figure 1, the resulting data set is dominated by the results of one study, 12 shown in grey. In an attempt to compensate for this, in the original analysis 1 the results of the other studies were given three times as much weight in the fitting. In Figure 2, 95% JCRs are shown for the four studies 12,[17][18][19] that provided sufficient data to estimate them, while the kp values at 25°C from the 4 remaining studies [13][14][15][16] are represented as open squares, shown at an arbitrary EA. Additionally, the 95% JCR 1 for the fit to the combined, weighted data is shown in red. It is apparent from Figure 2 Please do not adjust margins Please do not adjust margins individual JCRs, and is consistent with only one of the point estimates of kp at 25°C. This is a strong indication that there is significant interlaboratory variation in the data. In other words, we should not expect two laboratories to converge on the same parameter estimates, no matter how many times they repeat the experiments. Small systematic differences between laboratories in equipment, operator technique and raw materials result in observable differences in the parameter estimates. As such, the statistical model for the kp data should include a term for interlaboratory variation (eq. 5):

ln 𝑘 p,𝑖𝑗 = (ln 𝑘 p0 + 𝜀 𝑖 ) -

(𝐸 A +𝜂 𝑖 ) R ( 1 𝑇 - 1 𝑇 0 ) + 𝛿 𝑖𝑗 (5)
In eq 5, kp,ij represents the j th measurement of kp from the i th study. The random experimental error is represented by δij, while εi and ηi represent the error in the parameters kp0 and EA, respectively, associated with all measurements from the i th study. These errors are all assumed to be drawn from normal distributions with means of 0 and variances of V(δi), V(ε), and V(η). Note that V(δi) may vary from one study to another. The question then becomes: how can we estimate the Arrhenius parameters while taking into account the systematic interlaboratory variation? Equation 5suggests that the best estimates of V(ε) and V(η) will be obtained from the sum of squared differences between the parameter estimates of individual studies and the average of all studies. The parameter estimates of individual studies [12][13][14][15][16][17][18][19] of STY are grouped in Table 2, along with those of 8 additional studies [20][21][22][23][24][25][26][27] that were published after the original IUPAC paper, but which meet the conditions for inclusion. This gives a total of 16 independent estimates of ln kp at 25°C and 12 independent estimates of EA, with standard deviations of 0.082 and 0.77 kJ•mol -1 , respectively. The standard deviation in ln kp corresponds to a relative standard deviation (RSD) in kp of ±8%, which compares well with typical reproducibility standard deviations of up to 30% for molar masses obtained by SEC. 28,29 The mean of the individual studies provides an estimate of the Arrhenius parameters, which are found to be 4.46 ± 0.02 (ln kp) and (31.7 ± 0.2) kJ•mol -1 (EA), with a correlation coefficient, ρ, of 0.052 indicating that the uncertainties in the two values are essentially uncorrelated. The estimate for ln kp at 25°C is in good agreement with the original fit to the combined weighted data, but the estimate for EA is significantly (~3σ) lower. This is because less weight is given to the study with the most datapoints, 12 whose EA of 32.6 kJ•mol -1 dominated the combined fit presented in the original analysis. 1 As a result, the estimate for A is also significantly (~2σ) lower than the originally reported benchmark value, with a revised value of 10 7.49 ± 0.08 L•mol -1 •s -1 , compared to the originally published 1 value of 10 7.63 ± 0.06 L•mol -1 •s -1 .

In order to take the varying precision of the individual estimates into account, we then calculated the weighted average of the kp and EA estimates, weighting each pair of parameter estimates according to the sum of the estimated interlaboratory variance and the variance estimated from the individual study. For the studies reporting only one or two kp values, the reproducibility of the kp measurement was assumed to be ± 5%. We believe this to be a conservative estimate of precision in a single laboratory; in contrast to the larger variance in interlaboratory reproducibility of SEC, the intralaboratory reproducibility (i.e. the same sample reanalysed using the same equipment) is generally good. 28,29 Full details of the statistical treatment are given in the Supporting Information. This gave estimated values of 4.46 ± 0.03 for ln kp at 25°C (kp = (86 ± 3) L•mol -1 •s -1 ), and (31.7 ± 0.4) kJ•mol -1 for EA, values scarcely different to the unweighted average. The uncertainties in ln kp at 25°C and EA are only slightly correlated (ρ = 0.24). While the weighted average leads to a higher uncertainty in the estimate of EA, it remains significantly (~2σ) lower than the originally published value of 32.5 kJ•mol -1 . A 95% JCR 30,31 can be calculated for these parameters, but is significantly larger than the originally published JCR due to the severe reduction in degrees of freedom. The original dataset contains only 4 independent estimates of EA, compared to the 61 data points used to calculate the original JCR, while 12 estimates of EA are available using the expanded dataset. As a result, the critical value of the F-distribution used to obtain the JCR is taken from the F2,3 distribution (original dataset) or F2,11 distribution (expanded dataset) as 19.1 or 8.0, respectively, as opposed to the F2,59 distribution, which gives a critical value of 6.4. § Figure 3a shows the JCR for the expanded dataset (in blue), overlaid on the JCRs of the individual studies. This JCR encompasses significant portions of most of the individual JCRs and point estimates of kp. Figure 3b shows the comparison between the originally published JCR (red), the JCR recalculated from the original data of reference 1 taking into account the interlaboratory variation (orange), and the JCR obtained from the expanded dataset (blue). Gratifyingly, when interlaboratory variation is taken into account, the inclusion of additional studies results in a negligible change in the parameter estimates, while significantly improving their precision. While the revised parameter estimates fall outside the originally determined JCR; the original parameter estimates remain within the revised JCR, indicating that these original estimates remain consistent with the experimental data.

Extension to all benchmarked monomers

This approach to determining the JCR while taking into account interlaboratory variation works well when there are at least 5 independent studies. However, the majority of the IUPAC benchmark studies comprised 2 to 4 laboratories. With so few independent data points, the 95% JCRs become unfeasibly large. Thus, we sought an alternative method to estimate the typical interlaboratory variation in estimation of Arrhenius parameters for propagation rate coefficients by PLP. Assuming that the error was roughly constant regardless of the monomer being studied, we calculated the pooled interlaboratory covariance matrix of all 13 monomers for which Please do not adjust margins Please do not adjust margins benchmark Arrhenius parameters are available as the weighted average of the variances and covariances of the individual parameter estimates for each monomer, according to equation 7:

𝑠 𝑝𝑜𝑜𝑙𝑒𝑑 2 = ∑ (𝑁 𝑖 -1)𝑠 𝑖 2 𝑖 ∑ (𝑁 𝑖 -1) 𝑖 (7)
The pooled variance included 81 independent estimates of ln kp and 71 independent estimates of EA. This gave a standard deviation for the interlaboratory error in ln kp at 25°C of 0.08, and in EA of 1.4 kJ•mol -1 , with a correlation coefficient of 0.04. This can be used to calculate a 95% JCR for the interlaboratory error, which is shown in Figure 4. One study 32 on butyl acrylate was excluded from this calculation as its estimate of EA deviated substantially from the mean of the remaining studies of the same monomer. It should be noted that this study was carried out over a very small temperature range (5-25°C), which may explain the imprecision in the estimate of EA in this case. A second study 33 on methyl acrylate was also excluded as the results were obtained at high pressure and extrapolated back to ambient pressure. The deviations of the parameter estimates of each individual study from the mean EA and ln kp at 25°C for the appropriate monomer are also shown in Figure 4. Of the 71 points included, all but 3 fall within the estimated 95% JCR, in line with expectation (5% of 71 ≈ 4). Thus, we recommend that this estimate of uncertainty be applied to PLP studies from individual laboratories that have not yet been independently replicated. Applying this estimate of the interlaboratory variation to the case of STY, we arrive at essentially the same parameter estimates as before: ln kp at 25°C of 4.46 ± 0.02, EA of (31.8 ± 0.5) kJ•mol -1 , and a covariance of -0.37 kJ•mol -1 , corresponding to a correlation coefficient of -0.04. The standard error in ln kp is slightly reduced relative to the previous estimate, while that of EA is slightly greater. Full details of these calculations are given in the Supporting Information.

T (K) (b) (a) (c)
This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 7

Please do not adjust margins Please do not adjust margins parameter estimates are shown for comparison in Figure 5 for kp at 25°C and EA (Figure 5a) as well as A and EA (Figure 5b). For practical purposes there is little difference between the JCRs obtained for the Arrhenius parameters of STY from the pooled variance data for all monomers or from the styrene studies only. Both cover the same range of EA (30-33 kJ•mol -1 ) and a similar range of ln kp at 25°C (4.34-4.58 vs 4.37-4.54, equivalent to kp ranges of 77-98 and 79-94 L•mol -1 •s -1 , respectively). Both include the original parameter estimates of 4.44 and 32.5 kJ•mol -1 , although both new parameter estimates fall outside the originally published 95% JCR. Finally, the estimated uncertainty in kp tallies well with the intuition of experienced researchers in the field, who typically estimate an uncertainty of ± 10% in kp values obtained by PLP. The 95% confidence bands for kp as a function of temperature in the range of 260-365 K are shown in Figure 5c. If the pooled interlaboratory variance is used to calculate the confidence band, the uncertainty in kp is roughly ± 6% near 298 K, rising to approximately ± 10% at the extremities of the investigated temperature range. A slightly broader confidence band is obtained if only the styrene studies are used, with uncertainties of 8-15% depending on the temperature. Confidence bands for kp of all benchmarked monomers can be found in the Supporting Information. Similar calculations were carried out for all the monomers for which IUPAC benchmark data has been published, leading to the 95% JCRs and parameter estimates shown in Table 1 and Figure 6. For most monomers, the change in parameter estimates is minimal, however changes > 1σ in EA were obtained for styrene (a difference of 0.7 kJ•mol -1 ), and in kp at 25°C for butyl methacrylate, glycidyl methacrylate, butyl acrylate and methyl acrylate. The maximum difference between original and revised values was 7% (kp of GMA at 25°C) Comparing the JCRs shown in Figure 6, the trend towards higher kp for methacrylates with longer side chains is evident, while the EA seems independent of the side chain for all methacrylates except dodecyl methacrylate. Likewise the family-like behavior of the acrylates is clear, with butyl acrylate showing a higher kp than methyl acrylate but a similar activation energy.

Notes on data selection

Styrene. The data set from reference 1 (8 studies [12][13][14][15][16][17][18][19] ) was augmented with 8 additional studies [20][21][22][23][24][25][26][27] published between 1996 and 2006, and covering a temperature range of 18-120°C. Analysis was by SEC using polystyrene standards, with the exception of Willemse et al. 27 in which MALDI-TOF MS was used to obtain the MMDs.

Methyl methacrylate. The data set from reference 2 (7 studies 13,16,20,[34][35][36][37] ) was augmented with 12 additional studies 27,[38][39][40][41][42][43][44][45][46][47][48] published between 1997 and 2015, and covering a temperature range of -18-91.5°C. Analysis was by SEC using

Conclusions

Systematic differences between laboratories have a significant effect on the results obtained from PLP studies of the temperature dependence of the propagation rate coefficient of monomers in radical polymerization. These differences occurred despite significant efforts to standardize experimental protocols and must be taken into account when analyzing the results of interlaboratory studies. Failure to do so leads to significant underestimation of the uncertainty associated with the Arrhenius parameters, and may produce erroneous estimates when the experimental dataset is dominated by the results of a single study.

We have presented revised estimates of activation energies and pre-exponential factors (and their 95% JCRs) which explicitly account for the interlaboratory variation and also incorporate additional data sets published subsequently to the benchmark studies. In doing so, we have estimated the typical interlaboratory error as ± 0.08 for ln kp (equivalent to ± 8% in kp) and ± 1.4 kJ•mol -1 for EA, with a correlation coefficient of 0.04. This may be used as an estimate of the uncertainty of a single study of kp that has not been independently replicated.

The existence of systematic differences between laboratories underlines that in an interlaboratory study, it is better to have a relatively small number of results from many laboratories, rather than many results from a small number of laboratories.

Repetition of experiments can reduce the uncertainty in a single laboratory's result, but once this becomes small relative to the interlaboratory uncertainty, no further increase in precision is obtained from additional experiments. This applies equally to experiments where many data points are obtained from a single experiment (for example in the determination of reactivity ratios by tracking the change in monomer feed composition with conversion, or the determination of Mark-Houwink-Sakurada constants from online viscosimetry). In these cases, data points from a single experimental run should be assumed to be correlated, and multiple independent experiments should This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 9

Please do not adjust margins Please do not adjust margins be run in order to determine the uncertainty associated with the parameter estimates. In summary, we make the following recommendations:

1. The revised estimates of kp at 25°C, A, and EA and their associated uncertainties given in Table 1 replace the previously reported values [1][2][3][4][5][6][7][8] as IUPAC benchmarks. 2. The estimated interlaboratory error of ± 8% in kp and ± 1.4 kJ•mol -1 in EA should be assumed to apply to all PLP studies, and can provide a first estimate of the uncertainty in reported values when no independent replication is available. In addition, we note that the reported values for GMA show relatively poor agreement between replications, while the activation energy for dodecyl methacrylate is unusually low relative to the other methacrylates investigated. Further studies on these monomers would help to improve the accuracy of their parameter estimates.

Figure 1 .

 1 Figure 1. Dataset used to fit Arrhenius parameters for kp of styrene radical polymerization in bulk. Open circles (grey): data from reference 12; filled diamonds (grey): data from reference 17; filled circles (blue): data from reference 18; filled triangles (orange) data from reference 19; remaining points from references 13-16.

Figure 2 .

 2 Figure 2. 95% Joint confidence regions (JCR) for Arrhenius parameters of kp of styrene. In black, the JCRs corresponding to the individual studies of references 12 and 17-19. The open black squares represent 4 independent measures (references 13-16) of kp at 25°C, shown at an arbitrary EA. The red JCR corresponds to the original IUPAC benchmark fit (reference 1) to the combined, weighted data.

Figure 3 .

 3 Figure 3. 95% Joint confidence regions (JCRs) for individual studies (a) and pooled data (b) for EA and kp at 25°C of styrene radical polymerization in bulk. (a) JCRs corresponding to references 17-20 shown in grey, additional data shown in black. Open circles represent kp from single experiments (refs 13-16), shown at an arbitrary EA. Open square shows estimate of kp and EA from reference 20. As only two experiments were carried out, a JCR could not be determined. Blue JCR represents the JCR that corresponds to the average of all individual estimates. (b) Original IUPAC benchmark JCR and parameter estimates (red); revised JCR and parameter estimates based on original IUPAC dataset (orange); revised JCR and parameter estimates using extended dataset (blue).

6 |

 6 J. Name., 2012, 00, 1-3 This journal is © The Royal Society of Chemistry 20xx

Figure 4 .

 4 Figure 4. Interlaboratory variation in estimates of ln kp at 25°C and EA, and pooled 95% JCR. The point in red (reference 32) was considered an outlier and excluded from the study. Standard errors in ln kp at 25°C and EA are 0.08 and 1.4 kJ•mol -1 , respectively, with a correlation coefficient of 0.05

Figure 5 .

 5 Figure 5. (a) Comparison of parameter estimates and 95% JCRs for EA and ln kp at 25°C of radical polymerization of styrene in bulk obtained using estimates of interlaboratory variation obtained from studies on styrene (blue) or from pooled studies of all benchmarked monomers (green). Original IUPAC benchmark JCR and parameter estimates (red) are shown for comparison. (b) The same JCRs shown for the preexponential factor, A, and EA. (c) 95% confidence bands for kp of radical polymerization of styrene in bulk from 262 K to 364 K.

Figure

  Figure 6. Revised 95% JCRs for monomers for which IUPAC benchmark values are available. MMA: methyl methacrylate; EMA: ethyl methacrylate; BMA: butyl methacrylate; iBoMA: isobornyl methacrylate; DMA: dodecyl methacrylate; GMA: glycidyl methacrylate; cHMA: cyclohexyl methacrylate; BnMA: benzyl methacrylate; VAc: vinyl acetate; MAA: methacrylic acid; MA: methyl acrylate; BA: butyl acrylate.

8 |

 8 J. Name., 2012, 00, 1-3 This journal is © The Royal Society of Chemistry 20xx

Table 2 .

 2 Estimates of kp and EA from individual PLP studies of styrene radical polymerization in bulk, and statistics for the population of parameter estimates. Numbers in parentheses represent the standard error in the final digits: e.g. 4.344(50) represents 4.344 ± 0.050

	Reference	Ln kp at 25°C	EA (kJ•mol -1 )	n c
	13	4.344(50 b )		1
	14	4.431(50 b )		1
	15	4.364(50 b )		1
	16	4.522(50 b )		1
	12	4.497(07)	32.6(0.2)	45
	18	4.361(20)	31.1(0.9)	4
	19	4.538(11)	31.5(0.2)	4
	20 a	4.570(50 b )	31.3(1.7 b )	2
	21 a	4.399(39)	30.9(1.7)	10
	22 a	4.390(07)	31.8(0.3)	24
	23 a	4.482(06)	32.1(0.4)	4
	24 a	4.579(16)	33.1(1.1)	6
	25 a	4.372(19)	30.7(0.6)	4
	26 a	4.495(40)	31.9(0.6)	10
	27 a	4.547(09)	32.4(0.3)	3
	Variance	0.00666	0.585
	Covariance	0.00322
	Standard	0.082	0.76
	deviation
	Correlation	0.052
	coefficient
	(ρ) d
	Mean	4.46(2)	31.7(2)

a study published subsequent to original IUPAC STY paper. b assuming a standard error in kp of ± 5% c number of data points reported in each study d correlation coefficient, 𝜌 = 𝐶𝑜𝑣(ln 𝑘 𝑝 ,𝐸 𝐴 ) 𝜎 ln 𝑘𝑝 𝜎 𝐸 𝐴

  Please do not adjust marginsPlease do not adjust margins poly(MMA) standards, with the exception of Willemse et al.,27 in which MALDI-TOF MS was used to obtain the MMDs, and Gruendling et al.,47 in which coupled SEC/ESI-MS was used to determine accurate MMDs. Differences between kp values obtained by these methods and those obtained by SEC with MMA calibration were small compared to the typical interlaboratory variation. Analysis was by SEC using universal calibration, with the exception of Willemse et al,67 in which MALDI-TOF MS was used to obtain the MMDs. One early study32 was excluded from the analysis due to the unusually low reported activation energy. Methacrylic acid (15% in water). The data set from reference 5 (2 studies68,69 ) was reanalysed. No additional data were found that complied with the IUPAC guidelines for determination of kp using PLP. Methyl acrylate. The data set from reference 7 (5 studies33,[START_REF] Manders | [END_REF]67,[70][71] ) was reanalysed. Analysis was by SEC using universal calibration, with the exception of Willemse et al,67 in which MALDI-TOF MS was used to obtain the MMDs. The results of reference 33, extrapolated to ambient pressure from experiments at high pressure, were excluded. No additional data was found that complied with the IUPAC guidelines for determination of kp using PLP. Vinyl acetate. The data set from reference 8 (6 studies35,[START_REF] Manders | [END_REF][72][73][74][75] ) was reanalysed. No additional data were found that complied with the IUPAC guidelines for determination of kp using PLP.

	Ethyl methacrylate. The data set from reference 3 (3
	studies 44,49,50 ) was augmented with 1 additional study 47
	published in 2010, and covering a temperature range of 0-
	91.5°C. In this study, coupled SEC/ESI-MS was used to obtain
	the MMDs.
	Butyl methacrylate. The data set from reference 3 (4
	studies 44,49-51 ) was augmented with 4 additional studies 47,52-54
	published between 2004 and 2016, and covering a temperature
	range of 0-91.5 C. Analysis was by SEC using universal
	calibration, with the exception of Gruendling et al., 47 for which
	coupled SEC/ESI-MS was used to determine accurate MMDs.
	Differences between kp values obtained by this method and
	those obtained by SEC with universal calibration were small
	compared to the typical interlaboratory variation.
	Dodecyl methacrylate. The data set from reference 3 (3
	studies 49-51 ) was reanalysed. No additional data were found that
	complied with the IUPAC guidelines for determination of kp
	using PLP.
	Cyclohexyl methacrylate. The data set from reference 4 (3
	studies 55-57 ) was reanalysed. No additional data were found that
	complied with the IUPAC guidelines for determination of kp
	using PLP.
	Glycidyl methacrylate. The data set from reference 4 (2
	studies 55,57 ) was augmented with an additional study 58
	published in 2008. This study covered a temperature range of
	60-195°C, but only results from 60-90°C were added to the
	dataset in order to avoid contamination of the kp data from
	depropagation. In this study, SEC with both light scattering
	detection and universal calibration were used to obtain MMDs.
	The light scattering results were used in the current reanalysis.
	While the 3 studies gave consistent estimates of kp at 25°C, they
	differed quite significantly in their estimate of EA, ranging from
	20.3 to 26.2 kJ•mol -1 . As a result, the 95% JCR shown in Figure 6
	may underestimate the true uncertainty in EA for this monomer,
	and further studies would be helpful.
	Benzyl methacrylate. The data set from reference 4 (3
	studies 44,55,56 ) was augmented with 1 additional study 59
	published in 2011, and covering a temperature range of 14-
	72°C.
	Isobornyl methacrylate. The data set from reference 4 (2
	studies 55,56 ) was reanalysed. No additional data were found that
	complied with the IUPAC guidelines for determination of kp
	using PLP. It should be noted that these estimates were not
	considered a benchmark as the Mark-Houwink constants for
	iBoMA were not independently determined for each data set.
	Butyl acrylate. The data set from reference 6 (6 studies 32,60-64 )
	was augmented with 3 additional studies 65-67 published
	between 2008 and 2017, and covering a temperature range of -
	25-70°C.

Acknowledgements

This work was supported in part by the IUPAC Machine Readable Database Project. We acknowledge very helpful discussions with Alex van Herk.

Author Contributions

SH: Conceptualization, data curation, formal analysis, writingoriginal draft, writing -review and editing SB, RAH, TJ, GTR: conceptualization, writing -review and editing

Conflicts of interest

There are no conflicts to declare.

Notes and references