Ana Paula Amorim
email: apamorim@math.uminho.pt

Jacobo De Uña-Álvarez

Luís Meira-Machado

Presmoothing the transition probabilities in the illness-death model

Keywords: Kaplan-Meier, Markov condition, Multi-state models, semiparametric censorship

 estimator based on presmoothing. Consistency is established. We investigate the finite sample performance of the new estimator through simulations. Data from a study on colon cancer are used for illustration purposes.

Introduction

Multi-state models [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF][START_REF] Meira-Machado | Multi-state models for the analysis of time to event data[END_REF] are the most common models used for the description of longitudinal survival data. A multi-state model is a model for a stochastic process, which is characterized by a set of states and the possible transitions among them. [START_REF] De Uña-Álvarez | Strong consistency of presmoothed Kaplan-Meier integrals when covariables are present[END_REF][START_REF] Yuan | Semiparametric censorship model with covariates[END_REF] Iglesias-Pérez and de Uña-Álvarez, 2008). Recently, an application of presmoothing to the estimation of the bivariate distribution of censored gap times has been provided too (de [START_REF] De Uña-Álvarez | A semiparametric estimator of the bivariate distribution function for censored gap times[END_REF]. In this paper we will propose presmoothed estimators of the transition probabilities p ij (s, t) in the scope of the illnessdeath model.

In order to illustrate our estimators using real data, we consider data from one of the first successful trials of adjuvant chemotherapy for colon cancer, which is freely available from the R survival package. In this study, 929 patients affected by colon cancer underwent a potential curative surgery. Unfortunately, some of these patients had residual cancer, which lead to the recurrence of disease and death (in some cases). Therefore, we may consider the recurrence as an associated state of risk, and use the so-called illnessdeath model with states "alive and disease-free", "alive with recurrence" (local-regional or metastases) and "dead". See Section 2 for a more formal definition of the model.

The organization of the paper is as follows. Section 2 introduces the notation, the new estimators and the main results. The finite sample performance of the proposed estimator is investigated via simulations in Section 3. In Section 4 we analyze the colon cancer data with the proposed methods. Main conclusions are reported in Section 5. Further illustration, complete simulation results and technical proofs are deferred to the web-only Appendix.

The estimator: main result

In this paper we consider the illness-death model depicted in Figure 1. In this model, all the subjects are in State 1 ('healthy') at time t = 0. At some future time, they will arrive at State 3 ('dead'), which is absorbing. In the meanwhile they may visit State 2 ('diseased') at some time point; or not, passing directly to State 3 without visiting State 2. Note that this multi-state model is progressive [START_REF] Hougaard | Analysis of multivariate survival data[END_REF], in the sense that past states can not be revisited. For this model the set of states is S = {1, 2, 3}, and the transitions allowed are 1 → 2, 1 → 3, and 2 → 3. Given two time points s < t, there are in essence three different transition probabilities to estimate: p 11 (s, t), p 13 (s, t), and p 23 (s, t). The two other transition probabilities (p 12 (s, t) and p 22 (s, t)) are easily obtained from p 12 (s, t) = 1 -p 11 (s, t) -p 13 (s, t) and p 22 (s, t) = 1 -p 23 (s, t). More examples of multi-state models can be found in books by [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF] and [START_REF] Hougaard | Analysis of multivariate survival data[END_REF], or in papers by Putter et al. (2007) and Andersen and Perme (2008).

Despite its potential, multi-state modeling is not used by practitioners as frequently as other survival analysis techniques. It is our belief that lack of knowledge of available software and non-implementation of the new methodologies in userfriendly software are probably responsible for this neglect. One important contribution to this issue was given by the R/S-PLUS survival package. Thanks to this package, survival analysis is no longer limited to Kaplan-Meier curves and simple Cox models.

Indeed, this package enables users to implement the methods introduced by Therneau and Grambsch (2000) for modeling multi-state survival data. In R (R Development Core

Team 2008), multi-state regression can also be performed using the msm package Let T ij be the potential transition time from state i to state j. This means that a subject not visiting state 2 will reach the 'dead' state at time T 13 , while this time will be T 12 + T 23 if he/she passes through state 2 before. We denote by ρ = I(T 12 ≤ T 13 ) the indicator of visiting state 2 at some time. Let Z = T 12 ∧ T 13 be the sojourn time in state 1, and let T = Z + ρT 23 be the total survival time of the process (up to reaching the absorbing state). We denote the censoring variable by C which is assumed to be independent of the process; finally, we put Z = Z ∧ C and T = T ∧ C for the censored versions of Z and T , and ∆ 1 = I(Z ≤ C) and ∆ = I(T ≤ C) for the respective censoring indicators. With this notation, the transition probabilities are written as p 11 (s, t) = P (Z > t) P (Z > s) , p 13 (s, t) = P (s < Z, T ≤ t) P (Z > s) , p 23 (s, t) = P (Z ≤ s, s < T ≤ t) P (Z ≤ s < T ) .

All these quantities involve expectations of particular transformations of the pair (Z, T ), S (ϕ) = E [ϕ (Z, T )] say. Thus we now discuss how these expectations can be empirically approximated from the data

Z i , T i , ∆ 1i , ∆ i , ∆ 1i ρ i , 1 ≤ i ≤ n ,
which are assumed to form a random sample of the vector Z, T , ∆ 1 , ∆, ∆ 1 ρ . Note that p 11 (s, t) and the denominator of p 13 (s, t) only involve the Z variable, and that they can be estimated by the ordinary Kaplan-Meier estimator of the sojourn time distribution in state 1. However, the remaining quantities cannot be estimated so simply.

Let T 1:n ≤ • • • ≤ T n:n denote the ordered T i 's, and let W i be the Kaplan-Meier weight attached to T i:n when estimating the marginal distribution of T from ( T i , ∆ i )'s. That is,

W i = ∆ [i:n] n -i + 1 i-1 j=1 1 - ∆ [j:n] n -j + 1
where ∆ [i:n] is the i th concomitant of T i:n . Here, ties within the censored or within the uncensored times are ordered arbitrarily, and ties among the uncensored and censored times are treated as if the former precede the later.

In the uncensored case we have W i = n -1 for each i. In [START_REF] Meira-Machado | Nonparametric estimation of transition probabilities in non-Markov illnessdeath model[END_REF] the following estimator of S (ϕ) was proposed:

S n (ϕ) = n i=1 W i ϕ( Z [i:n] , T i:n ).
where Z [i:n] is the concomitant of T i:n . Consider now the presmoothed version of S n (ϕ) given by

S n (ϕ; m n ) = n i=1 W i (m n )ϕ( Z [i:n] , T i:n ) where W i (m n ) = m n Z [i:n] , T i:n n -i + 1 i-1 j=1   1 - m n Z [j:n] , T j:n n -j + 1  
and where m n (z, t) stands for an estimator of the binary regression function

m(z, t) = P ∆ = 1| Z = z, T = t .
Since Z, T , ∆ are observable, the function m(z, t) can be estimated by standard methods. However, the naive construction of a smooth estimator for m(z, t) will generally fail. This is because the function m(z, t) will typically be discontinuous along the line t = z, that is, for those covariate values Z, T corresponding to individuals who are censored while being in state 1 or who suffer a direct transition 1 → 3 to the absorbing state.

In order to see this, note that for z < t we have

m(z, t) = P ∆ = 1| Z = z, T = t, ∆ 1 ρ = 1 ≡ m 1 (z, t),
that is, m 1 Z, T is the conditional probability of censoring on T given Z, T and given that transition 1 → 2 is observed (∆ 1 ρ = 1). However, when z = t we get

m(t, t) = P ∆ 1 = 1| Z = t, ∆ 1 ρ = 0 ≡ m 2 (t),
which is the conditional probability of observing 1 → 3 given Z = t (or T = t) and given that transition 1 → 2 is never observed. We implicitly assume that the events Z = T and ∆ 1 ρ = 0 are the same. This is reasonable unless there is a significative proportion of individuals with zero sojourn time in state 2. These formulae show that the functions m 1 and m 2 represent different binary regression problems and that they are based on disjoint subpopulations (according to the value of ∆ 1 ρ). Furthermore, the limit of m 1 (z, t) as z approaches to t does not coincide with m 2 (t) in reality. Figure 2 displays these functions for the colon cancer data, when estimated separately by two logistic models. The noise around m 1 (z, t) comes from the fact that the variable z is omitted from the plot while it is present in the model (although without reaching statistical significance, p-value=0.285). Both functions are clearly separated.

In summary, in order to construct m n (z, t) we propose that one estimate the functions m 1 (z, t) and m 2 (t) independently by fitting some smooth models, m 1n (z, t) and m 2n (t) say, so we finally have

m n (z, t) = m 1n (z, t) I(z < t) + m 2n (t)I(z = t), or m n Z i , T i = m 1n Z i , T i I( Z i < T i ) + m 2n ( Z i )I( Z i = T i ) = m 1n Z i , T i ∆ 1i ρ i + m 2n ( Z i )(1 -∆ 1i ρ i ).
The estimator m 1n (z, t) is based on the subsample {i : ∆ 1i ρ i = 1}, while m 2n (t) is computed from {i : ∆ 1i ρ i = 0}. The only condition we assume on these two functions is that they should approximate well their targets in a uniform sense; more specifically, set

U 1 : sup z,t |m 1n (z, t) -m 1 (z, t)| → 0
w. p. 1, q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q 0 500 1000 Since m(z, t) = m 1 (z, t) I(z < t) + m 2 (t)I(z = t), under U 1 and U 2 we have

U : sup z,t |m n (z, t) -m(z, t)| → 0 w. p. 1.
and hence Theorem 2.1 in de Uña-Álvarez and Rodríguez-Campos ( 2004) can be applied with some adaptation to the present context (see the Appendix). Conditions under which U 1 and U 2 hold are investigated in a number of papers, including Devroye (1978aDevroye ( , 1978b)), [START_REF] Mack | Weak and strong uniform consistency of kernel regression estimates[END_REF], and [START_REF] Härdle | Uniform consistency of a class of regression function estimators[END_REF]. Now we state our main result and the corresponding corollaries. Let H be the distribution function of T and let

τ H = inf {t : H(t) = 1}.
Theorem 1. Assume that H is continuous, that U 1 and U 2 hold, and that

E |ϕ (Z, T )| I (T ≤ τ H ) m(Z, T )(1 -H(T )) ρ < ∞ A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT is satisfied for some ρ > 0. Then, S n (ϕ; m n ) → S τ (ϕ) with probability 1, where S τ (ϕ) = E [ϕ(Z, T )I(T ≤ τ H )].
Theorem 1 is a proper adaptation of the Strong Law in Dikta (2000) to our scenario. We note that the result is not restricted to parametric presmoothing; the only thing one should have in mind is that condition U must be verified by the chosen estimator m n (z, t). Note also that, in general, one can not ensure that S τ (ϕ) and S (ϕ) will coincide; indeed, as always with censored data, one should not expect consistency beyond the upper bound of the censoring distribution, because there is no sampling information regarding the lifetime there. As a particular case, we have S τ (ϕ) = S (ϕ) if the support of T is contained in that of C.

The proof to Theorem 1 is similar to that of Theorem 2.1 in de Uña-Álvarez and Rodríguez-Campos ( 2004); here, the role of their covariate vector is played by the sojourn time in State 1 (Z), while the total time T up to reaching the absorbing State 3 is taken as the 'response'. However, since Z is also subject to right-censoring, the results in the referred paper (see also [START_REF] Stute | Consistent estimation under random censorship when covariables are present[END_REF][START_REF] Stute | Consistent estimation under random censorship when covariables are present[END_REF][START_REF] Stute | The strong law under random censorship[END_REF] do not directly apply here. Note that, since C is assumed to be independent of (Z, T ), the identifiability conditions H1 and H2 in de Uña-Álvarez and Rodríguez-Campos (2004) automatically hold. In our setup, these conditions read H1. T and C are independent H2. P (T ≤ C|Z, T ) = P (T ≤ C|T ) which clearly follow from the independence between the censoring time and the process. See the web-appendix for details. Now, we come back to our initial goal of estimating the transition probabilities p ij (s, t). Recall that p 11 (s, t) can be estimated by the ordinary Kaplan-Meier based on the Z i , ∆ 1i 's. In order to introduce some presmoothing, we recommend to replace the ∆ 1i 's by some smooth fit to the binary regression function P ∆ 1 = 1| Z = z (see e.g. [START_REF] Dikta | On semiparametric random censorship models[END_REF]. Now, we focus on the estimation of p 13 (s, t) and p 23 (s, t). Write

p 13 (s, t) = P (s < Z, T ≤ t) P (s < Z) = E [ϕ s,t (Z, T )] P (s < Z) ,
where ϕ s,t (u, v) = I (u > s, v ≤ t). Introduce the presmoothed estimator

p 13 (s, t) = S n (ϕ s,t ; m n ) P (s < Z)
where P (s < Z) stands for a consistent estimator (e.g. Kaplan-Meier) of P (s < Z).

Similarly, we have

p 23 (s, t) = P (Z ≤ s, s < T ≤ t) P (Z ≤ s < T ) = E [ ϕ s,t (Z, T )] E [ϕ s (Z, T )]
where

ϕ s,t (u, v) = I (u ≤ s, s < v ≤ t) and ϕ s (u, v) = I (u ≤ s < v).
Therefore, in this case, we estimate the transition probability through

p 23 (s, t) = S n ( ϕ s,t ; m n ) S n (ϕ s ; m n ) .
We have the following Corollary.

Corollary 1. Assume that the conditions in Theorem 1 hold for the special ϕ-functions ϕ s,t , ϕ s,t and ϕ s . Then, for any consistent estimator P (Z > s) of P (Z > s) we have with probability 1 p 13 (s, t) → p τ 13 (s, t) and p 23 (s, t) → p τ 23 (s, t), where p τ 13 (s, t) = P (T ≤ t, T ≤ τ H / Z > s) and p τ 23 (s, t) = P (T ≤ t / Z ≤ s < T, T ≤ τ H ).

Corollary 1 is an immediate consequence of Theorem 1. As for the Theorem, consistency can not be ensured in general. However, when the support of C contains that of T we have p τ 13 (s, t) = p 13 (s, t) and p τ 23 (s, t) = p 23 (s, t). In particular this may happens whenever τ H = ∞.

Remark In practice, when n is small, it may happen p 13 (s, t) > 1 and/or p 11 (s, t) + p 13 (s, t) > 1. When any of these inequalities occurs, we propose the modification p 13 (s, t) = 1 -p 11 (s, t), which ensures p 13 (s, t) ≤ 1 and p 11 (s, t) + p 13 (s, t) ≤ 1. With this remark in mind, we always have p 12 (s, t) = 1 -p 11 (s, t) -p 13 (s, t) ≥ 0. For moderate or large sample sizes this problem disappears.

Simulation study

In this Section we investigate the performance of the proposed estimators p ij (s, t) through simulations. More specifically, the estimators p 11 (s, t), p 13 (s, t) and p 23 (s, t) introduced in Section 2 are considered.

To simulate the data in the illness-death model, we separately consider the subjects passing through State 2 at some time (that is, those cases with ρ = 1), and those who directly go to the absorbing State 3 (ρ = 0). For the first subgroup of individuals (ρ = 1), the successive gap times (Z, T -Z) are simulated according to the bivariate distribution

F 12 (x, y) = F 1 (x)F 2 (y) [1 + θ {1 -F 1 (x)} {1 -F 2 (y)}]
where the marginal distribution functions F 1 and F 2 are exponential with rate parameter 1. This corresponds to the so-called Farlie-Gumbel-Morgenstern copula, where the single parameter θ controls for the amount of dependency between the gap times. The parameter θ was set to 0 for simulating independent gap times, and also to 1, corresponding to 0.25 correlation between Z and T -Z. This simulated scenario is the same as that described in [START_REF] Lin | Nonparametric estimation of the gap time distributions for serial events with censored data[END_REF] and de Uña-Álvarez and [START_REF] De Uña-Álvarez | A simple estimator of the bivariate distribution function for censored gap times[END_REF]. For the second subgroup of individuals (ρ = 0), the value of Z is simulated according to an exponential with rate parameter 1. In summary, the simulation procedure is as follows:

Step 1. Draw ρ ∼ Ber(p) where p is the proportion of subjects passing through State 2.

Step 2. If ρ = 1 then:

(2.1) V 1 ∼ U (0, 1) , V 2 ∼ U (0, 1) are independently generated;

(2.2)

U 1 = V 1 , A = θ (2U 1 -1)-1, B = (1 -θ (2U 1 -1)) 2 +4θV 2 (2U 1 -1) (2.3) U 2 = 2V 2 / √ B -A (2.4) Z = ln (1/ (1 -U 1 )) , T = ln (1/ (1 -U 2 )) + Z If ρ = 0 then Z = ln (1/ (1 -U (0, 1))).
Situations with p = 1 corresponds to the three-state progressive model, in which a direct transition 1 → 3 is not allowed. In our simulation we consider p = 0.7. An independent uniform censoring time C is generated, according to models U [0, 4] and U [0, 3]. The first model results in 24% of censoring on the first gap time Z, and in 47% of censoring on the second gap time T -Z,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

for those individuals with ρ = 1. The second model increases these censoring levels to 32% and about 57%, respectively.

After some algebra, it is seen that the function

m 1 (z, t) = P ∆ = 1| Z = z, T = t, ∆ 1 ρ = 1 is written as m 1 (z, t) = 1 1 + η 1 (z, t) , where η 1 (z, t) = λ G (t) λ 1 T |Z=z (t|z)
and where λ G (.) and λ 1 T |Z=z (.|z) stand respectively for the hazard rate of the censoring variable and the hazard rate of T given Z = z under restriction

ρ = 1. Note that λ G (t) = 1/(τ G -t) when C ∼ U τ G ] and that λ 1 T |Z=z (t|z) is given by λ 1 T |Z=z (t|z) = 2 + 4 exp(-t) -2 exp(-z) -2 exp(-t + z) 2 + 2 exp(-t) -2 exp(-z) -exp(-t + z) if θ = 1,
being 1 when θ = 0. The function m 1 (z, t) belongs to the logistic family with some preliminary transformation of the conditioning variables. To be more specific we have (for β 0 = 0 and β 1 = 1)

m 1 (z, t; β) = 1 1 + exp(β 0 + β 1 ln(η 1 (z, t)))
.

This is the parametric model we fit to m 1 (z, t) in the simulations. The β parameter in model m 1 (.; β) is estimated via maximization of the conditional likelihood of the ∆ i 's given the ( Z i , T i )'s, for those subjects with ∆ 1 ρ = 1 (see e.g. [START_REF] Dikta | On semiparametric random censorship models[END_REF][START_REF] Dikta | The strong law under semiparametric random censorship models[END_REF]. The same estimation criterium is used for the other presmoothing functions (m 0 and m 2 ) in this section. For m 2 (t

) = P ∆ 1 = 1| Z = t, ∆ 1 ρ = 0 , we have m 2 (t) = 1 1 + η 2 (t) , where η 2 (t) = λ G (t) λ 0 Z (t)
and where λ 0 Z (t) stands for the sub-hazard function of Z restricted to ρ = 0, namely λ 0 scenario with θ = 1 (dependent transition times). The row labeled m corresponds to presmoothing with the true function m(z, t) = P (∆ = 1| Z = z, T = t), which is unrealistic because this function will generally be unknown. However, this row represents a 'gold standard' the other methods can be compared to. The row labeled with m(.; β, γ) corresponds to a semiparametric estimator which is obtained using a presmoothing based on a parametric family which contains the true m. Specifically, we consider a logistic model with the preliminary transformation of the conditioning variables Z = z, T = t shown before. Results for p 11 (s, t) and p 13 (s, t) are shown in the web-appendix. Similarly, for p 11 (s, t) and for the denominator of p 13 (s, t) we also perform logistic presmoothing for the function m 0 (z) = P ∆ 1 = 1| Z = z , with the variable Z transformed by -ln(τ G -Z) (so the parametric family contains the true m 0 (z)).

Z (t) = P (Z = t, ρ = 0|Z ≥ t) = 1 -p.

Statistic

In order to investigate the robustness of the proposed estimator with respect to miss-specifications of the binary regression family, we considered also presmoothing via standard logistic models, without any preliminary transformation of the transition times. This is labeled with m(., ξ) in Tables 2 and3. Note that the true m and the true m 0 do not belong to this parametric family. Finally, we also report the results pertaining to the estimator in [START_REF] Meira-Machado | Nonparametric estimation of transition probabilities in non-Markov illnessdeath model[END_REF], which corresponds to the situation with no presmoothing at all. This is labeled in the Tables as KM.

Some expected features are clearly seen in Tables 2 and3. For example, we see that the (integrated) MSE, bias and variance of p23 (s, t) decrease with an increasing sample size, while they increase with the censoring degree. The best performance is attained by the estimator which makes use of the true m, which was expected. However, in practice one must estimate the function m. The lowest errors among the realistic versions of the estimators correspond to the estimator based on the correctly specified parametric family, m(.; β, γ). Finally, we see that the presmoothed estimator based on the wrong parametric model m(.; ξ) is still (much) better than KM; the practical message is that it is worthwhile doing some presmoothing even when we are not completely sure about the parametric family.

Compared to the estimator without any presmoothing (KM), it is proven that the relative efficiency of the estimators based on presmoothing is always above 1. In special cases, the relative deficiency of the Kaplan-Meier estimator is below 50%; this occurs for larger values of s, where the censoring effects are stronger. This supports the belief that the relative benefits of presmoothing will be seen more clearly in the presence of large censoring degrees.

Although we restrict the integrated bias, variance and mean square error (MSE) to the interval [s, 3], we verified that, in both settings, the enlargement of this interval favors the estimator based on presmoothing (results not shown). This happens because higher levels of censoring are expected in the right tail of the distribution. We also run simulations for mutually independent transition times T ij . The results (shown in the web-appendix) also revealed advantages on the use of presmoothing.

Colon cancer study

For illustration, we apply the proposed methods of Section 2 to data from a large clinical trial on Duke's stage III patients, affected by colon cancer, that underwent a curative surgery for colo-rectal cancer [START_REF] Moertel | Levamisole and fluorouracil for adjuvant therapy of resected colon carcinoma[END_REF]. In this study, from the total of 929 patients, 468 developed recurrence and among these 414 died. 38 patients died without recurrence. The rest of the patients (423) remained alive and disease-free up to the end of the followup. As mentioned in the Introduction recurrence can be expressed as an intermediate event which can be modeled using an illness-death model.

Using the Cox proportional hazards model, we verified that the transition rate from state 2 to state 3 is affected by the time spent in the previous state. This allowed us to conclude that the Markov assumption may be unsatisfactory for the colon cancer data set. In this section we will present estimated transition probabilities calculated using the new approach, based on presmoothed Kaplan-Meier weights and the estimators of [START_REF] Meira-Machado | Nonparametric estimation of transition probabilities in non-Markov illnessdeath model[END_REF]. Neither one of the approaches assume the process as being Markovian.

In Figure 3 we illustrate differences between the estimated transition probabilities p ij (s, t), 1 ≤ i ≤ j ≤ 3 based on presmoothing the Kaplan-Meier weights and the estimator corresponding to no presmoothing (KM; [START_REF] Meira-Machado | Nonparametric estimation of transition probabilities in non-Markov illnessdeath model[END_REF]. The presmoothed estimator was obtained by standard logistic regression for both m 1 and m 2 . The value s was chosen to be the 75th percentile of the sojourn time in state 1 (s = 1549 days). From this figure we conclude that the new estimator have more jump points (corresponding to patients with censored values of the total time) but with smaller steps. The number of jump points and the size of the steps are related to the censoring degree and to the sample size. We can also verify that both methods provide similar point estimates for small time values. Departures between both estimated curves can be more appreciated for larger time values where the censoring effects are stronger. In summary, the new approach provides more reliable curves with less variability, specially at the right tail of the lifetime distribution. Other values of s reported similar results.

Conclusions and final remarks

In this paper we have introduced new estimators for the transition probabilities of a censored illness-death model. The new estimators are based on a preliminary estimation (presmoothing) of the probability of censoring for the total time, given the available information. This idea has been used before in univariate survival analysis, but its application to multi-state survival data is complicated since new problems arise. More explicitly, the illness-death model involves presmoothing functions which are discontinuous, so a naive estimation approach fails.

We have derived the consistency of the proposed estimators. The consistency result is not restricted to parametric presmoothing, but it also includes the possibility of using some nonparametric estimators to this end. The finite sample performance of the introduced estimators was investigated through simulations. The main conclusion is that presmoothing leads to improved estimators, even when there is some miss-specication in the parametric family assumed for the presmoothing function. The relative benefits of presmoothing are more clearly seen in the heavily censored case. The new method has been illustrated using data from a colon cancer study.

The new estimators for the transition probabilities are consistent regardless the Markov condition. This is interesting because real problems are often far from markovianity and therefore the consistency of the time-honored Aalen-Johansen estimator can not be ensured. To this regard, one may think about the methods introduced here as a remarkable improvement (in the

Figure 2 :

 2 Figure 2: Illness-death model.

(

  continuous-time Markov and hidden Markov multi-state models), the changeLOS package(Wrangler et al. 2006) implements the Aalen-Johansen estimator for general multi-state models, and the etm package has recently enabled the transition matrix to be computed, along with a covariance estimator.This paper describes the R-based p3state.msm package's capabilities for analyzing survival data from an illness-death model. It extends existing semi-parametric regression capabilities included in many statistical software programs, such as R, S-PLUS, SAS, etc. Moreover, p3state.msm enables several quantities of interest to be estimated, such as transition probabilities, bivariate distribution function, etc. In

Figure 1 :

 1 Figure 1: Illness-death model: the three states (boxes) and the possible transition among them (arrows).

Figure 2 :

 2 Figure 2: Presmoothing functions m 1 and m 2 estimated by logistic models vs. T (variable Z not shown). Colon cancer data.

and U 2

 2 : sup t |m 2n (t) -m 2 (t)| → 0 w. p. 1.

Figure 3 :

 3 Figure 3: Estimated transition probabilities for p ij (s, t) with s = 1549 based on the Kaplan-Meier weights (dashed line) and based on presmoothed Kaplan-Meier weights (solid line). Colon cancer data.

Table 1 :

 1 Summary statistics measuring bias, variance and mean square error.

		Definition	Estimator
	Integrated absolute bias Integrated variance Integrated MSE	t 1 s |bias(s, t)|dt t 1 s var(p ij (s, t))dt t 1 s M SE(p ij (s, t))dt	t 1 t=s | bias(s, t)|δ t 1 t=s var(p ij (s, t))δ t 1 t=s M SE(p ij (s, t))δ

Acknowledgement. The authors tank the AE and a referee for suggestions which have improved the presentation of the paper. Work supported by the Grants MTM2008-03129 of the Spanish Ministerio de Ciencia e Innovación and 10PXIB300068PR of the Xunta de Galicia. Support from the INBIOMED project of the Xunta de Galicia (DXPCTSUG, Ref. 2009/063) is also acknowledged. Luis F. Meira-Machado acknowledges financial support by Grant PTDC/MAT/104879/2008 (FEDER support included) of the ACCEPTED MANUSCRIPT Portuguese Ministry of Science, Technology and Higher Education. Ana P. Amorim and Luis Meira-Machado acknowledge financial support provided by the research Centre of Mathematics of the University of Minho through the FCT Pluriannual Funding Program.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Similarly as above, we fit the logistic model m 2 (t; γ) = 1 1 + exp(γ 0 + γ 1 ln(η 2 (t)))

to estimate the function m 2 (t) in the simulations. As before, this logistic model has the true presmoothing function m 2 as a special case (γ 0 = 0,γ 1 = 1).

The aim of this simulation study is to compare the estimator by [START_REF] Meira-Machado | Nonparametric estimation of transition probabilities in non-Markov illnessdeath model[END_REF] and the new estimator based on presmoothing ideas. In order to measure the estimates' relative performance, we computed the integrated absolute bias, integrated variance and the integrated Mean Square Error (MSE) of the estimates. For each simulated setting we derived the analytic expression of p 11 (s, t), p 13 (s, t) and p 23 (s, t) so that the bias and the MSE of the estimator could be examined. Sample sizes 50, 100 and 200 were considered. In each simulation, K = 1000 samples were generated.

Let pk ij (s, t) denote the estimated transition probability based on the kth generated data set. For each fixed (s, t) we obtained the mean for all generated data sets, pij (s, t)

We then computed the pointwise estimates of the bias, variance and MSE as:

To summarize the results we also calculated the integrated absolute bias, integrated variance and the integrated MSE, defined in Table 1. We fixed the values of s using the quantiles 0.25, 0.5 and 0.75 of the exponential distribution with rate 1. The results obtained in Table 2 and 3 were obtained by numerical integration on the interval [s, t 1 ] with t 1 = 3, taking a grid of step δ = 0.05.

In Tables 2 and3 we report the results for the integrated absolute bias, integrated variance and the integrated MSE attained by the proposed estimators for p 23 (s, t) when based on several presmoothing functions, in the


Geert Molenberghs

Geert Verbeke

Samuel Iddi

Pseudo-likelihood methodology for partitioned large and complex samples

Keywords: Asymptotic relative efficiency, Compound-symmetry, Small-sample relative efficiency

Large data sets, either coming from a large number of independent replications, or because of hierarchies in the data with large numbers of within-unit replication, may pose challenges to the data analyst up to the point of making conventional inferential methods, such as maximum likelihood, prohibitive. Based on general pseudo-likelihood concepts, we propose a method to partition such a set of data, analyze each partition member, and properly combine the inferences into a single one. It is shown that the method is fully efficient for independent partitions, while with dependent sub-samples efficiency is sometimes but not always equal to one. It is argued that, for important realistic settings, efficiency is often very high. Illustrative examples enhance insight in the method's operation, while real-data analysis underscores its power for practice.

Introduction

Contemporary statistical analysis is confronted with data sets of ever increasing dimensions, not only because the number of independent units may become very large, but also owing to hierarchies in the data with large numbers of within-unit replication. While computational and data-analytic resources have recently progressed tremendously, sheer size may be a limiting factor. Our method allows for data analysis that otherwise would be prohibitive, either due to extremely large sample size or a large amount of dependent replication.

Using pseudo-likelihood technology, we propose a method whereby a dataset can be chopped into portions, each of which is then analyzed separately, followed by combining inferences into a single set.

The method is inspired by work of Fieuws and Verbeke (2006) and Fieuws et al (2006), reviewed in Molenberghs and Verbeke (2005, Ch. 24), who constructed a tool for the analysis of high-dimensional longitudinal data.

Our focus is on proper partitions, where every data point is a member of one and exactly one of the, say M , subsamples. Our method could be generalized to situations with overlapping subsamples as well. However, proper partitions allow us to use maximum likelihood on each subsample, thus simplifying the combination of M inferences into a single one considerably.

A particularly interesting special case is when partitioning is done in terms of independent units (i.e., dependent data are not spread over different members subsample). In that case, the method is fully efficient in the sense that the Cramèr-Rao lower bound is reached. This is true, even though the actual estimator may differ from the one that would result were the data analyzed in full. These results are illustrated using univariate normal and univariate Bernoulli samples.

When correlated data are partitioned by splitting a sequence of dependent data across different subsamples, then full efficiency is sometimes, but not always, reached. An example of this setting is when for each subject in a longitudinal study, M n repeated measurements are taken, which are subsequently split into M sequences of length n. This results in M separate longitudinal sets of data, each one having shorter length. This situation is useful when the length of a sequence is prohibiting efficient computation. Unlike in the previous case, a dependent partition results. This case is illustrated using the simple, insightful example of a multivariate normal sample with constant mean and compound-symmetry covariance structure. It will be shown that such dependent partitioning can be fully efficient for some but not necessarily for all parameters. The latter remarks begs the questions as to the efficiency of the method. The advantage of our illustrative examples is that the asymptotic relative efficiency, and its small-sample counterparts, can be studied without difficulty.

Of course, the illustrative examples have limited use in practice. The method is devised for settings that genuinely poses computational challenges that can be alleviated by partitioning the data set.

We apply our method to a data comprising long binary repeated-measures sequences.

The paper is organized as follows. The motivating data are introduced in Section 2; their analysis is relegated to Section 6. The methodology for our method is presented in Section 3, with attention for general pseudo-likelihood concepts, the method of Fieuws and Verbeke, and our partitioned-sample case, with emphasis on both independent and dependent partitions. The aforementioned illustrations are described in Section 4 and efficiency scrutinized in Section 5.

Motivating Case Studies

The motivating case study encompasses data collected by the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA; [START_REF] Mathei | Evidence for a substantial role of sharing of injection paraphernalia other than syringes/needles to the spread of hepatitis C among injecting drug users[END_REF]. It consist of annual serological surveys, providing information about the hepatitis C virus (HCV) and human immunodeficiency virus (HIV) status and related risk factors from the 20 Italian regions in the period 1998-2006. Respondents are drug users who sought help in the drug treatment centers. In each year, large numbers of drug users are tested for their HCV and HIV statuses. The resulting dataset is a sequence of binary outcomes for each of the 20 Italian regions over the years. Not all regions are present at all years. We will focus on the HCV data. The maximum number of respondents is 15,401 (average 3866.61), and the maximum of positive tests is 10,875 (average 2578.12). The aim is to investigate the change in HCV over time, i.e., whether year-effects are present in the profiles. Table 3 in the Supplementary Material lists the overall HCV prevalences for each of the 20 regions. Figure 1 in the Supplementary Material displays the profiles for each of the 20 regions.

Methodology

Standard Pseudo-likelihood

Using [START_REF] Arnold | Pseudolikelihood estimation: some examples[END_REF], we introduce pseudo-likelihood, the principal idea of which is to replace a numerically challenging joint density by a simpler function assembled from suitable factors.

Consider a sample of size N with repeated measures sequences of length n. Let S be the set of all 2 n -1 vectors of length n consisting solely of zeros and ones, with each vector having at least one non-zero entry. Denote by y (s) i the subvector of y i corresponding to the components of s that are non-zero. The associated joint density is f s (y (s) i ; β). To define a pseudo-likelihood function, one chooses a set δ = {δ s |s ∈ S} of real numbers, with at least one non-zero component. This is a general definition and the precise choice depends on inferential goals and the specific form of the model being considered. A common choice is pairwise likelihood where δ s = 1 for sequences s with exactly two components equal to 1, and zero otherwise. Other choices include the set of full conditionals, ensuring every outcome in a sequence is conditioned upon all other outcomes in the sequence of repeated measures. Details can be found in Molenberghs and Verbeke (2005, Ch. 9, 12, 21, 22, 24, and 25). The log of the pseudo-likelihood is then

p = N i=1 s∈S δ s ln f s (y (s) i ; β).
(1)

Adequate regularity conditions have to be invoked to ensure that (1) can be maximized by solving the pseudo-likelihood (score) equations, the latter obtained by differentiating the logarithmic pseudolikelihood and equating its derivative to zero. These regularity conditions are spelled out in the Supplementary Material. In particular, when the components in (1) result from a combination of marginal and conditional distributions of the original distribution, then a valid pseudo-likelihood function results. In particular, the classical log-likelihood function is found by setting δ s = 1 if s is the vector consisting solely of ones, and 0 otherwise. Broadly, a pseudo-likelihood is valid if composed of marginal and conditional densities, derived from the full density describing the entire sequence of measurements, thereby allowing for repetition and weighting. More details can be found in [START_REF] Varin | On composite marginal likelihoods[END_REF], [START_REF] Lindsay | Composite likelihood methods[END_REF], and [START_REF] Joe | On weighting of bivariate margins in pairwise likelihood[END_REF]. Note that [START_REF] Joe | On weighting of bivariate margins in pairwise likelihood[END_REF] use weighting for reasons of efficiency in pairwise likelihood, similar in spirit to [START_REF] Geys | A note on the comparison of pseudo-likelihood and generalized estimating equations for marginal odds ratio models[END_REF], but differently from its use here, which focuses on bias correction when data are incomplete. Another important reference is [START_REF] Cox | A note on pseudolikelihood constructed from marginal densities[END_REF].

Let β 0 the true parameter. Under suitable regularity conditions [START_REF] Arnold | Pseudolikelihood estimation: some examples[END_REF][START_REF] Geys | Pseudo-likelihood modelling of multivariate outcomes in developmental toxicology[END_REF][START_REF] Aerts | Topics in Modelling of Clustered Data[END_REF], it can be shown that maximizing (1) gives a consistent, asymptotically normal estimator β N so that √ N ( βNβ 0 ) converges in distribution to

N p [0, I 0 (β 0 ) -1 I 1 (β 0 )I 0 (β 0 ) -1 ]. ( 2
)
Precise statements and additional discussion are given in Appendix A of the Supplementary Material section.

3.2 Pseudo-likelihood for Partitioned Samples

Background

Fieuws and Verbeke ( 2006) and Fieuws et al. (2006) proposed a pseudo-likelihood-based method to fit mixed models to high-dimensional longitudinal data. Their method is reviewed in Molenberghs and Verbeke (2005, pp. 470ff). Precisely, when a large number of longitudinal sequences are modeled simultaneously, standard (restricted) maximum likelihood becomes prohibitive. As an alternative, they propose fitting corresponding mixed models to each pair of outcomes. Hence, if there are M longitudinal sequences per subject, M (M -1)/2 pairs ensue. The difference with standard pseudolikelihood, as reviewed in Section 3.1, is that the models are fitted to each pair separately, whereas (1) assembles all contributions into a single pseudo-likelihood function, ensuring that all parameters are estimated only once. In the pairwise approach, the mean parameters governing, for example, the first longitudinal sequence, are estimated M -1 times, because the first second is paired with the second, the third, and so on, up to the M th. Nevertheless, these authors are able to cast their method in the general pseudo-likelihood context. To see this, they first observe that fitting all bivariate models is equivalent to maximizing the function

p (θ) ≡ p (y 1i , y 2i , . . ., y M i |θ) = r<s (y ir , y is |θ rs ), (3)
where y mi is sequence m = 1, . . . , M for subject i = 1, . . . , N , θ is the overall parameter vector and θ rs is the parameter vector pertaining to pair (r, s). To proceed, one temporarily ignores that some of the vectors θ rs have common elements, i.e., assuming that all vectors θ rs are completely distinct. In (3), θ results from stacking all M (M -1)/2 pair-specific parameter vectors θ rs . The actual parameter vector of interest is θ * , the set of non-redundant parameters is θ.

Evidently, (3), is of the form (1) and hence their pairwise fitting procedure fits within the general framework of pseudo-likelihood. Conveniently, the set of parameters in θ rs is treated pair-specific, which allows separate maximization of each term in the pseudo log-likelihood function (3).

Because the pairwise approach fits within the pseudo-likelihood framework, an asymptotic multivariate normal distribution for θ can be derived, using the general pseudo-likelihood theory presented in Section 3.1:

√ N( θ -θ) approx.
∼ N (0, I -1 0 I 1 I -1 0 ). To pass from θ to θ * , Fieuws and Verbeke (2006) take averages of all available estimates for that specific parameter, implying that θ * = A θ for an appropriate linear combination matrix A. Hence, inference for the elements in θ * will be based on

√ N( θ * -θ * ) = √ N(A θ -A θ) approx. ∼ N (0, A I -1 0 I 1 I -1 0 A). (4)
We will make use of these ideas to apply pseudo-likelihood estimation to partitioned samples, partioned pseudo-likelihood (PPL) for short.

Independent Subsamples

Let us modify the ideas, reviewed in Section 3.2.1, to the case where a given sample, deemed too large, is broken into m = 1, . . . , M subsamples, each of size n, such that N = M • n. While it is possible to let the subsample size vary (with then N = M m=1 n m ), to avoid cluttering notation we focus on the equal subsample size case. The likelihood for sample m takes the form

p m (θ m ) = n i=1 (y mi |θ m ), (5)
where (•) refers to the likelihood one would consider were the mth subsample the entire set of data.

Further, y mi is the ith subject in sample m. Note that all θ m are equal to θ * , the parameter vector of interest. Note that θ in Section 3.2.1 now takes the form (θ * , θ * , . . . , θ * ). The overall estimator is naturally defined as:

θ * = 1 M M m=1 θ m . ( 6
)
Owing to the independence of the subsamples, the off-diagonal blocks in the matrices I 0 and I 1 are all zero. Because each of the M contributions to (5) pertains to a genuine likelihood, the diagonal blocks in I 0 and I 1 are identical up to the sign, because

H θ ≡ I 0 (θ m , θ m ) = E ∂ 2 m (θ m ) ∂θ T m ∂θ m = -E ∂ m (θ m ) ∂θ m T . ∂ m (θ m ) ∂θ m = -I 1 (θ m , θ m ).
Using the above results, the asymptotic variance (4) for (6) takes the form

var( θ * ) = 1 M H -1 θ . (7)
While ( 7) uses the expected information, one can also employ the observed information matrices, denoted by H θ,m and derived as the second derivative of the pseudo-likelihood function, for each of the subsamples:

1 M 2 M m=1 H -1 θ,m . ( 8
)
These derivations are based upon using

A = 1 M (I, . . ., I) (9)
in ( 4), where I is an identity matrix with dimensions equal to the length of the vector θ * .

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Dependent Subsamples

In the above, the overall sample was partitioned into M independent parts, consisting each of n subjects. A different, important area of application is when, say longitudinal or otherwise hierarchical samples would be partitioned by splitting each outcome vector Y i into M sub-vectors:

(Y 1i , . . ., Y M i ) , each consisting of n measurements. Three remarks are in order. First, as before, it is perfectly possible to let n vary with subsample, and even from subject to subject, but for ease of exposition we will keep it constant throughout this section. Second, partitioning can be generalized to that of the previous and the current section combined. Third, the difference with Section 3.2.1 is that we do not make multiple use of the same portions within an outcome sequence; should this be desirable, this aspect can then be brought into the picture without trouble. The methodology in the dependent case is entirely similar to that of the dependent case with the sole but important exception that the off-diagonal blocks in I 1 are generally non-zero. Indeed, the (m, m ) block of I 1 takes the form

I 1 (θ m , θ m ) = E ∂ m (θ m ) ∂θ m T . ∂ m (θ m ) ∂θ m .
As a result, the straightforward expressions ( 7) and ( 8) no longer apply. Rather, general expression (4) should be used, with A as in (9).

Illustrations

In this section, we provide some simple but insightful illustrations in the form of a univariate normal sample, Bernoulli sampling, and a compound-symmetry multivariate normal sample. While a genuine data analysis is relegated to Section 6, these three simple cases have the advantage of establishing three qualitatively distinct situations: (1) No difference between likelihood (ML) and PPL. In the normal sample as well as in the Bernoulli sample, parameterized using a probability, there is no difference between either. The same is true for the mean parameter in the compound-symmetry case;

(2) Different estimator, same precision. When the Bernoulli experiment is characterized by the logit instead, a different estimator emerges, with nevertheless the same asymptotic variance; Different estimator, precision loss. ML and PPL differ for both variance components (measurement error, random-intercept variance) in the compound-symmetry case, in terms of both the estimators as well as their asymptotic variance. Each of these three cases will be sketched next, with details relegated to the Supplementary Material, where appropriate.

A Univariate Normal Sample

Assume Y mi ∼ N (µ, σ 2 ). The PPL takes the form:

p (µ 1 , . . ., µ M , σ 2 1 , . . . , σ 2 M ) ∝ M m=1 n i=1 - 1 2 ln(σ m ) - 1 2 y mi -µ m σ m 2 . ( 10
)
Recall that all µ m and all σ m are equal across m, but in the process of deriving the estimator, each subsample is given its own parameter vector. Focusing on the mean parameter, we obtain µ m = y mi .

The mean estimator becomes

µ = 1 M M m=1 1 n n i=1 y mi = µ, (11)
where µ refers to the PPL estimator, as opposed to the ML estimator µ. The second derivative of ( 10) is -n/σ 2 m . Given equality of the parameters, we can write I 0 (µ m , µ m ) = -I 1 (µ m , µ m ) = nσ 2 and, with now A = 1 M (1, . . ., 1) , we obtain var( µ) = σ 2 /(M n) = var( µ). Based upon similar calculations, the asymptotic variance for σ 2 becomes: var( σ 2 ) = 2σ 4 /(M n) = var( σ 2 ). Note that we have treated both estimators independent rather than joint, because µ and σ 2 are independent, a property carrying over onto both I 0 and I 1 , so that also in the PPL case these are effectively two scalar-parameter problems.

Hence, as stated earlier, for the univariate normal case, the PPL method produces the same estimators and the same asymptotic variances as the classical one obtained with ML. Given the linearity of all expressions involved, this result is also true for the observed information. When the sample sizes would vary across sub-samples, then (11) has to be weighted appropriate to maintain equality.

A Univariate Bernoulli Sample

Assume that Y mi is Bernoulli with parameter π or, equivalently, logit α = ln[π/(1 -π)]. The PPL assumes the form:

p (π 1 , . . . , π M ) = M m=1 [z m ln π m + (n -z m ) ln(1 -π m )] , (12)
with z m = n i=1 y mi , the number of successes in the mth subsample, for which π m is the corresponding copy of the success probability π. It then follows that π m = z m /n and hence

π = 1 M n M m=1 z m = π,
where π refers to the ML estimator and π is the PPL counterpart. Further

E ∂ 2 m ∂π 2 m = - n π m (1 -π m ) , leading to var( π) = [π(1 -π)]/(M n) = var( π).
Again, in analogy with the univariate normal case above, there is no difference between the PPL and ML cases.

Switching to the logit parameterization, ( 12) is replaced with

p (α 1 , . . ., α M ) = M m=1 [z m α m -n ln (1 + e αm )] , (13)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

leading to the asymptotic variances

α m = ln z m n -z m , and E ∂ 2 m ∂α 2 m = -n e αm (1 + e αm ) 2 ,
and hence, the variance for the logit becomes:

var( α) = 1 M 2 n M m=1 1 π m (1 -π m ) . ( 14
)
In case we use the expected value, we obtain:

var( α) = 1 M n 1 π(1 -π) ,
exactly the same as what would be obtained if the sample were not partitioned but rather conventional ML were used. This is interesting, because the point estimator

α = 1 M M m=1 α m is different from the ML version α = ln[z/(M n -z)],
where z is the overall number of successes out of a total of M n Bernoulli trials.

Multivariate Normal Compound Symmetry

Consider the compound symmetry model

Y i = (Y 1i , . . . , Y M i ) ∼ N (1µ, σ 2 I + dJ), ( 15
)
where 1 is a vector of ones of appropriate length, I is the identity matrix, and J is a matrix of ones.

The dimensions of 1, I, and J will be suppressed from notation when no confusion can arise. We will present the ML and PPL point and precision estimators. Details are relegated to Appendix B in the Supplementary Material.

To simplify notation, let us consider the pseudo-likelihood contribution for a given arbitrary subsample, m, which comes down to the log-likelihood for a sample of N subjects with sequences of length n. Recall that the full sequence length is M n and that the sequences have been split into M equal portions. The log-likelihood takes the form:

(µ m , σ 2 m , d m ) = - 1 2 N i=1 ln σ 2n + nσ 2(n-1) d m +(Y mi -µ m 1) 1 σ 2 m I - d m σ 2 m + nd m (Y mi -µ m 1) . ( 16
)
As in the univariate normal case, the mean µ m on the one hand and the variance components σ 2 m and d m on the other can be treated separately, owing to their functional and statistical independence.

Considering µ m first and setting the first derivative of ( 16) equal to zero leads to the conventional:

µ m = 1 N n N i=1 n j=1 Y ij . (17)
From ( 17) we immediately deduce:

µ = 1 N nm M m=1 N i=1 n j=1 Y mij . ( 18
)
which is the overall mean. Evidently, this is also the estimator obtained, should one analyze the sample as a whole, i.e., µ = µ. Similarly, turning to the variance components, it follows that

σ 2 m = 1 N n(n -1) n N i=1 Z mi Z mi - N i=1 Z mi J n Z mi , (19)
d m = 1 N n(n -1) N i=1 Z mi J n Z mi - N i=1 Z mi Z mi , (20)
where Z mi = (Y mi -µ1 n ) and J n is an n × n matrix of ones. From ( 19) and ( 20) it follows that:

σ 2 = 1 M N n(n -1) n M m=1 N i=1 Z mi Z mi - M m=1 N i=1 Z mi J n Z mi , (21)
d = 1 M N n(n -1) M m=1 N i=1 Z mi J n Z mi - M m=1 N i=1 Z mi Z mi . (22)
Evaluating ( 19) and ( 20) for the case that the sample is analyzed at once, but then evidently for sequences of length M n, leads to:

σ 2 = 1 M N n(M n -1) M n N i=1 Z i Z i - N i=1 Z i J M n Z i , (23)
d = 1 M N n(M n -1) N i=1 Z i J M n Z i - N i=1 Z m Z m , (24)
where

Z i = (Y i -µ1 M n ) and J M n is an M n × M n matrix of ones. Thus, σ 2 = σ 2 and d = d.
Consider precision estimation. For this, the (block-diagonal) matrix I 0 and the (non-block-diagonal)

I 1 need to be derived. Details are found in Supplementary Material B. As stated earlier, the mean can be treated separately from the variance components. Turning to µ, given that µ = µ, it is not surprising that both asymptotic variances are equal:

var( µ) = var( µ) = σ 2 + M nd M N m . ( 25
)
Focusing on the variance components, we find:

var( σ 2 ) = 2σ 4 M N n -M N , (26)
var( σ 2 ) = 2σ 4 M N n -N , (27) var
( d) = 2 M N n σ 4 M n -1 + 2dσ 2 + M nd 2 , (28)
var( d) = 2 M N n σ 4 n -1 + 2dσ 2 + M nd 2 . ( 29
)
5 Efficiency

Asymptotic Relative Efficiency

In Section 4, three distinct illustrations were used. In the univariate normal and univariate Bernoulli cases (probability and logit scale), and for the mean parameter in the compound-symmetry case, the asymptotic relative efficiency (ARE), defined as the variance ratio of the ML over the PPL estimator, when the overall sample size tends to infinity, equals one. For the logit case, this may not be the case for small samples, a point to which we will in Section 5.2. This is no longer the case for the variance components in the compound-symmetry case. Indeed, for σ 2 , we find:

ARE(σ 2 ) = M n -M M n -1 . ( 30
)
Three remarks are in place. First, if M = 1, there is no partitioning, and evidently ARE reduces to 1. Second, if n = 1, the ARE ánd the variances are problematic. In particular, if n = 1 and M > 1, the ML variance ( 27) is well defined, whereas PPL variance ( 26) is undefined, because for the PPL each of the subsamples is univariate, from which σ 2 and d cannot be disentangled. Third, if n approaches infinity, with M bounded, the ARE approaches one. The above remarks also hold for the random-intercepts variance d, even though the ARE is slightly more tedious:

ARE(d) = 1 M n-1 + 2 d σ 2 + M n d σ 2 2 1 n-1 + 2 d σ 2 + M n d σ 2 2 = (1-ρ) 2 M n-1 + 2ρ(1 -ρ) + M nρ 2 (1-ρ) 2 n-1 + 2ρ(1 -ρ) + M nρ 2 , ( 31
)
where ρ = d/(σ 2 + d), the intraclass correlation. We will study these expressions in what follows.

The ARE for σ 2 has a simple structure, and nicely shows the efficiency loss related to increasing the number of partitions of the original sample. When the sequence length per component, n, is large relative to the number of sub-samples, efficiency loss will be modest. For example, with M = 5 pieces and n = 20, the ARE is around 95%. The sequence length is realistic because, for very sort sequences, the need to revert to sub-divided samples is usually not present.

ARE expression (31), shows that the two key quantities are M , the number of subsamples in the partition, on the one hand, and the intraclass correlation ρ on the other. Evidently, ARE(d) approaches 1 when M approaches 1; for large M , there is a 'partitioning penalty.' When ρ is small, the penalty is large, because then (31) approaches ARE(d|ρ = 0) = (n -1)/(M n -1). For example, when the correlation is close to 0 and M = 5, estimating the random-intercept variance is effectuated at a This is a trivial extension of the methodology presented so far. For convenience, we ensured that portions were roughly equal, because exact equality could not be achieved. Now, weighting matrix A to take the form (9) may not be fully optimal; although this does not affect the validity of the method, efficiency may be affected slightly. Second, the said efficiency in the dependent case is more affected when parameterization (32) is used, for the parameters α 1 -α 8 , than with parameterization (33), with efficiency largely unaffected. Also, under (32), the overall intercept parameter α 0 is unaffected by partitioning. Note that parameters β 1 -β 9 have the meaning of an intercept parameter as well, be it for a specific year, whereas α 1 -α 8 are contrasts between a given year and the last one. Hence, in some parameterizations, the efficiency loss may be spread out relatively evenly across parameters, whereas others may segregate parameters that are largely unaffected from others that are more severely impacted. Even for strongly affected parameters, α 1 -α 8 , it seems to be more a consequence of partitioning as such, rather than of the number of components, M , in the partition. This is not general; counterexamples are given by ( 31) and ( 30), which are smooth, continuous

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Table 1 :

 1 HCV data. Parameter estimates and standard errors for independent partitioning. ForIn line with Section 3.2.2, results are virtually identical across values of M . In the final analysis, dependent samples are created by sub-dividing the sequences into M = 1, 2, 5, 10, 15 parts. Results are in Table2. Several observations can be made. First, sub-samples are not all of equal size.

	M = 2, 4, 15 the proper empirically corrected standard errors are followed by their (inappropriate)
	purely model-based counterparts.
	Par. M = 1 (ML)	M = 2	M = 4
		Parameterization (32)
	α 0	0.592(0.112) 0.598(0.111) 0.593(0.108)
	α 1	0.223(0.011) 0.213(0.011) 0.243(0.012)
	α 2	0.209(0.011) 0.202(0.011) 0.215(0.011)
	α 3	0.288(0.011) 0.287(0.011) 0.300(0.012)
	α 4	0.179(0.011) 0.175(0.011) 0.170(0.011)
	α 5	0.106(0.011) 0.099(0.011) 0.095(0.011)
	α 6	0.114(0.011) 0.104(0.011) 0.106(0.011)
	α 7	0.072(0.011) 0.062(0.011) 0.068(0.011)
	α 8 -0.037(0.011) -0.043(0.011) -0.049(0.011)
	σ	0.501(0.079) 0.493(0.078) 0.459(0.076)
		Parameterization (33)
	β 1	0.815(0.113) 0.811(0.111) 0.836(0.108)
	β 2	0.801(0.113) 0.800(0.111) 0.808(0.108)
	β 3	0.880(0.113) 0.886(0.111) 0.894(0.108)
	β 4	0.771(0.113) 0.773(0.111) 0.763(0.108)
	β 5	0.698(0.112) 0.697(0.111) 0.689(0.108)
	β 6	0.706(0.112) 0.702(0.111) 0.699(0.108)
	β 7	0.664(0.112) 0.660(0.111) 0.662(0.108)
	β 8	0.555(0.113) 0.556(0.111) 0.544(0.108)
	β 9	0.592(0.112) 0.598(0.111) 0.593(0.108)
	σ	0.501(0.079) 0.493(0.078) 0.459(0.076)

Table 2 :

 2 HCV data. Parameter estimates and standard errors for dependent partitioning. For M =

	1, 2, 5, 10, 15, the proper empirically corrected standard errors are followed by their (inappropriate)
	purely model-based counterparts.
	Par. M = 1 (ML)	M = 2	M = 5	M = 10	M = 15
				Parameterization (32)
	α 0	0.592(0.112) 0.592(0.119;0.080) 0.592(0.119) 0.593(0.119) 0.593(0.119;0.030)
	α 1	0.223(0.011) 0.223(0.077;0.011) 0.223(0.077) 0.223(0.077) 0.223(0.077;0.011)
	α 2	0.209(0.011) 0.209(0.070;0.011) 0.209(0.070) 0.209(0.070) 0.209(0.070;0.011)
	α 3	0.288(0.011) 0.288(0.063;0.011) 0.288(0.063) 0.288(0.063) 0.288(0.063;0.011)
	α 4	0.179(0.011) 0.179(0.061;0.011) 0.179(0.061) 0.179(0.061) 0.179(0.061;0.011)
	α 5	0.106(0.011) 0.106(0.055;0.011) 0.106(0.055) 0.106(0.055) 0.106(0.055;0.011)
	α 6	0.114(0.011) 0.114(0.051;0.011) 0.114(0.051) 0.114(0.051) 0.114(0.051;0.011)
	α 7	0.072(0.011) 0.072(0.054;0.011) 0.072(0.054) 0.072(0.054) 0.072(0.054;0.011)
	α 8 -0.037(0.011) -0.037(0.033;0.011) -0.037(0.033) -0.037(0.033) -0.037(0.033;0.011)
	σ	0.501(0.079) 0.501(0.079;0.056) 0.500(0.079) 0.498(0.079) 0.496(0.079;0.021)
				Parameterization (33)
	β 1	0.815(0.113) 0.815(0.096;0.080) 0.815(0.096) 0.815(0.096) 0.816(0.096;0.030)
	β 2	0.801(0.113) 0.801(0.116;0.080) 0.802(0.116) 0.802(0.116) 0.802(0.116;0.030)
	β 3	0.880(0.113) 0.880(0.127;0.080) 0.880(0.127) 0.881(0.127) 0.881(0.127;0.030)
	β 4	0.771(0.113) 0.771(0.112;0.080) 0.772(0.112) 0.771(0.113) 0.771(0.113;0.030)
	β 5	0.698(0.112) 0.698(0.121;0.080) 0.699(0.121) 0.699(0.121) 0.699(0.121;0.030)
	β 6	0.706(0.112) 0.706(0.119;0.080) 0.706(0.119) 0.707(0.119) 0.707(0.119;0.030)
	β 7	0.664(0.112) 0.664(0.131;0.080) 0.665(0.131) 0.666(0.131) 0.666(0.131;0.030)
	β 8	0.555(0.113) 0.555(0.118;0.080) 0.556(0.118) 0.556(0.118) 0.557(0.118;0.030)
	β 9	0.592(0.112) 0.592(0.119;0.080) 0.593(0.119) 0.593(0.119) 0.593(0.119;0.030)
	σ	0.501(0.079) 0.501(0.079;0.079) 0.500(0.079) 0.498(0.079) 0.496(0.079;0.021)

Acknowledgment

We acknowledge financial support from IAP research Network P6/03 of the Belgian Government (Belgian Science Policy).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

mere 20% of the ML efficiency. This is qualitatively among the worst scenarios. Of course, it actually corresponds to the situation where the component d is unimportant in the first place. Luckily, on the other hand, the penalty vanishes for d large, relative to σ 2 , because ARE(d|ρ = 1) = 1.

Small-sample Relative Efficiency for the Logit Parameter

To assess the small-sample impact when there is asymptotic equivalence, but not equality of the estimators, we conduct limited simulations for a logit estimated from univariate Bernoulli samples, as in Section 4.2. We compare partition-based variance ( 14) with its single-sample counterpart, when the true logit is 0 (probability of 0.5). The number of partition components is M ∈ {1, 5, 10, 20, 50} and the sample size per component is n ∈ {20, 50, 100, 500}. Results are in Table 4 of the Supplementary Material. We see that the SSRE is very high. Given that the ARE has been shown to equal 1, it is not a surprise that the SSRE is very close to 1 as well when n is large relative to M . Only when the number of partition components is large relative to the sample size per partition component, does the SSRE shrink somewhat, but the worst case remains at a comfortably high 90%.

Analysis of Case Study

We analyze the data of Section 2. For the ith region in year j, let Z ij be the number of independent reported cases of HCV out of n ij . Assume a binomial model with success probability:

with the indicator for year defined as T ij = 1 if the year during which the measurement is taken equals k = 1, . . ., 9, and 0 otherwise. The study spans 9 years, with the last year as reference.

The model is of the generalized linear mixed model type [START_REF] Molenberghs | Models for Discrete Longitudinal Data[END_REF], with a region-specific effect b i ∼ N (0, σ 2 ). To study the effect of parameterization on efficiency, consider the equivalent formulation:

For each of the parameterizations, we consider three approaches.

In the first analysis, the data are analyzed as if measurements are independent (i.e., omitting the random effect). Results are summarized in Table 5 of the Supplementary Material. This assumption is unrealistic; the analysis is included to make a few points about the methodology only. It confirm the results for the (independent) Bernoulli case. The difference with Section 4.2 is that here we assume logistic regression models. Not unexpectedly, the parameter estimates and standard errors are the same, up to four decimal places for all of the values of M = 1, 2, 5, 10, 15. As a second analysis, the data are split into M independent subsamples, for M = 1, 2, 4. Results are in Table 1.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

functions of M . Arguably, efficiency will be more at risk for parameters that have a within-subject meaning. This agrees with the fact that, in the compound-symmetry case (4.3), µ is unaffected by partitioning, whereas the variance components, especially d, can be affected considerably stronger.

Third, it is unwise to replace the empirically corrected standard errors, based on I -1 0 I 1 I -1 0 with their model-based counterparts, using I -1 0 only, as can be seen when comparing the appropriate empirically corrected standard errors with their model-based counterparts, in Tables 1 and2. This is in line with other uses of pseudo-likelihood [START_REF] Molenberghs | Models for Discrete Longitudinal Data[END_REF].

Concluding Remarks

We have presented a convenient, simple, pseudo-likelihood based method to partition large data sets to facilitate estimation. Partitions into dependent as well as independent subsamples have been studied. In the independent subsample case, full efficiency can be reached. This is not always so for dependent samples, but in important realistic settings high to very high efficiency can be obtained.

Especially in such cases, it is recommendable to ensure the sub-sample size is not too small. Needless to say that large sets of data are extremely common in current-day empirical research, including but not limited to large surveys, microarray experiments, consumer databases, etc. When confronted with choice, independent sub-samples are preferred, because of efficiency. Dependent sub-samples become indispensable whenever analyzing the sequences in full is practically prohibitive.

Finally, whether splitting in sub-sampling is done and, if so, in how many components, strongly depends on what is practically and numerically feasible. Indeed, given the potential loss in efficiency, fewer components is better than more.

Apart from illustrative examples, real-data analysis has been employed to underscore the practical power of the method. While the illustrative examples allow for closed-form derivation of point and precision estimators, this is relevant from theoretical and illustrative perspectives only. In particular, it is is clear that efficiency loss in the dependent case is relative to a particular parameter, which is why the mean, variance, and correlation parameters all exhibit different behavior. Therefore, the data analyst is advised to study efficiency loss, perhaps using a simulation study designed after the real application at hand. Of course, there are cases where maximum likelihood is not even feasible, underscoring our method's use.

For data analysis, such as that based on the generalized linear mixed model in Section 6, closed-form derivations are neither feasible (because even conventional ML defies closed forms) nor necessary (because our methodology is presented in generic terms and can be applied whenever score-vector contributions are available, either in analytical or numerical form). SAS code used for the purpose of this article are available from the authors' web site.