A Nested Copula Duration Model for Competing Risks with Multiple Spells
Lo, Simon M. S. ; Mammen, Enno; Wilke, Ralf
Document Version
Accepted author manuscript
Published in:
Computational Statistics & Data Analysis
DOI:
10.1016/j.csda.2020.106986
Publication date:
2020
License CC BY-NC-ND
Citation for published version (APA):
Lo, S. M. S., Mammen, E., & Wilke, R. (2020). A Nested Copula Duration Model for Competing Risks with Multiple Spells. Computational Statistics & Data Analysis, 150, [106986].
https://doi.org/10.1016/j.csda.2020.106986 Link to publication in CBS Research Portal
General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.
Take down policy
If you believe that this document breaches copyright please contact us (research.lib@cbs.dk) providing details, and we will remove access to the work immediately and investigate your claim.
Download date: 06. Nov. 2022
A Nested Copula Duration Model for Competing Risks with Multiple Spells
Simon M.S. Lo
∗Enno Mammen
†Ralf A. Wilke
‡April 2020
Abstract
A copula graphic estimator for the competing risks duration model with multiple spells is presented. By adopting a nested copula structure the dependencies between risks and spells are modelled separately. This breaks up an implicit restriction of popular duration models such as multivariate mixed proportional hazards. It is shown that the dependence structure between spells is identifiable and can be estimated, in contrast to the dependence structure between competing risks. Thus, by allowing these two components to differ, the model is not identifiable. This is an important finding related to the general identifiability of competing risks models. Various features of the model are investigated by simulations and its practicality is illustrated by an application to unemployment duration data.
Keywords: Nested Archimedean copula, Multiple occurrences, Frailty
1 Introduction
Duration models are important for the analysis of data in many disciplines, including biostatis- tics, mechanical engineering and social sciences. As a running example we consider the duration
∗City University Hong-Kong, Department of Economics and Finance, E–mail: losimonms@yahoo.com.hk
†Heidelberg University, Institute for Applied Mathematics, E–mail: mammen@math.uni-heidelberg.de
‡Copenhagen Business School, Department of Economics, E–mail: rw.eco@cbs.dk
of unemployment in a model with competing risks and multiple spells. Competing risks mean unemployment is terminated by different exit types such as starting a job (risk 1) or withdrawal from the labour market (risk 2). Multiple spells or repeated occurrences correspond to having more than one unemployment period per unemployed person in the data. We suggest a new copula model that permits for flexible modelling of dependencies between competing risks and multiple spells.
In a competing risks model, observables are not the latent competing durations of different exit types, but only the minimum of them and the cause of risk. In our running example, if an unemployed was observed to withdraw from the labour market, the latent unemploy- ment duration terminated by finding a new job is incidentally censored and is unobserved.
As a consequence, the joint distribution of the risk-specific latent durations as well as the respective marginal distributions are not identifiable (Cox, 1962; Tsiatis, 1975). Identifiable are only cause-specific sub-distributions, which are the distributions of the minimum of the latent durations and the exit type. In our running example, these are the distributions of the unemployment durations that are terminated by risk 1 and risk 2 respectively.
Copulas are increasingly popular in statistics to model the joint distribution of multiple out- comes (e.g. Kole, et al., 2007; Joe, 2015; Oh and Patton, 2015; Klein and Kneib, 2016; Klein et al., 2019). The advantage of using a copula function for modeling the joint distribution is that it separates the dependence structure from the marginal distributions, so that the copula function can be modelled separately from the marginals. Zheng and Klein (1995) show that the latent marginals in the competing risks model can be nonparametrically identified by establishing the relationship between the latent marginals and the cause-specific sub-distributions provided that the copula function is known. It is called the copula graphic estimator, which is compatbile with a large variety of copula classes, including Archimedean, hierarchical Archimedean, vine copula, and factor copula models. For a review of different types of copula functions, see Joe (2015). The wide range of copula models provides flexibilities in modelling high-dimensional multivariate joint distributions in terms of the number of parameters, features of tail depen- dence (e.g. different degrees of dependence between extreme events), and different degrees of asymmetries (e.g. the dependence between two particular risks can be different from the dependence between another pair of risks). The Archimedean copula is mainly used for the copula graphic estimator, though, as it permits for flexible dependence structures, while also
taking into account the ease of implementation. Specifically, Rivest and Wells (2001) obtain a closed form relationship between the latent marginals and the cause-specific sub-distributions if the copula belonged to the Archimedean class. Similar linkages between the cause specific sub-distributions and the latent marginal distributions have been established in terms of hazard functions for Archimedean copula (Emura et al., 2019). Despite its computational tractability, the Archimedean class is rich enough to model different types of tail asymmetries, e.g. the Clay- ton copula exhibits lower tail dependence, the Gumbel copula exhibits upper tail dependence, while the Joe-Clayton copula allows for different degrees of lower and upper tail dependence.
At the same time it is parsimonious enough to be summarised by one to two parameters.
Popular alternatives to the copula graphic estimator include the assumption of independent risks (Edin, 1989; Mealli and Pudney, 1996; Lindstrom and Lauster, 2001; Burda et al., 2015;
Braun et al. 2010) and an assumed multivariate normal joint distribution (Malevergne and Sornette, 2003; D’ Addio and Rosholm, 2005; van den Goorbergh et al., 2005; Pennington- Cross 2010; MacKenzie and Spears, 2014; Smith, and Vahey, 2016; Agnello et al., 2019). These alternatives and the copula model share the common limitation that the choice of dependence structure is typically not guided by existing theoretical research. In our illustrative exam- ple regarding unemployment duration particularly, we are not aware of any economic theory that provides prior hypotheses on the choice of copula. Nevertheless, as a result of the non- identifiability aspect of the competing risks model, a certain degree of assumptions on the dependence structure must be made and these assumptions cannot be tested for. In this re- gard, the copula graphic estimator provides more flexibility than the alternatives as it allows the researcher to work with different copulas to fit the data. If the copula function were un- known, one could also conduct sensitivity analysis for the choices of the copula function. For instance, Lo and Wilke (2014) perform a series of Monte Carlo simulation studies and observe that estimation results are more sensitive to the choices of copula parameter(s) and less variant with respect to the choices of the Archimedean copula function.
The use of multiple spells data to address the non-identifiably of the competing risks model has also gained popularity in practice as the existence of several observations from one unit contributes additional information to the model. Honor´e (1993) shows that the mixed propor- tional hazard model (MPHM), as a specific class of duration models, is identified under much weaker restrictions in presence of multiple spells than if there were only one spell per unit.
This finding has since then triggered an extensive economic literature for multiple spells data (for example Horowitz and Lee, 2004; Kalwij, 2010; Arranz and Canto, 2012; Carrasco and Garcia-Perez, 2015). The key component to model dependencies between risks and between spells in the MPHM is a random frailty term. Frailty corresponds to some omitted variables that play a role for the latent durations. In our running example, the frailty term could be the degree of risk-aversion of an unemployed person, which is usually unobserved. More risk- adverse individuals tend to have shorter unemployment duration due to a sooner start of a new job in the trained occupation or in a different profession. More risk-adverse individuals also tend to have faster job-taking times for different unemployment spells. In this regard, the unobserved risk aversion determines the (positive) dependency between risks (types of jobs) as well as between spells. It is therefore an implicit assumption in the MPHM that the risks and the spells share the same dependence structure. Even though it is difficult to establish a 1-1 link between the mixed proportional hazard model and the copula based duration model, Lo et al. (2017) establish a linkage between the frailty term in the MPHM and the copula function by showing that the single spell MPHM (Hausmann and Han, 1990) is a special case of the copula model. Other contributions that consider copula duration models and frailty are Ha et al. (2019) and Wang et al. (2020). We reason in this paper that the frailty in multiple-spell MPH model is also related to a multiple-spell variant of the copula model.
We contribute to this literature by considering a nested copula model that accounts for (possibly different) dependencies between the competing risks and between multiple spells.
This relaxes the implicit assumption of the MPHM that the risks and the spells share the same dependency structure (Honor´e, 1993, Horny and Picchio, 2010), which can be important in some applications. In the context of our running example, unobservable frailty could also be the motivation of the unemployed. It is intuitive that the unobserved motivation causes a positive dependency between spells, as a more motivated unemployed person tends to have shorter unemployment duration for all her unemployment spells. However, whether the unobserved motivation causes a positive or a negative dependency between risks is unclear. For instance, higher levels of motivation may speed up time to employment but may slow down the time to drop out of the labour force. This means, motivation may have opposite effects on the risks, which translate into positive or negative risk dependency. So, it is more difficult to anticipate the direction of risk dependency. And more importantly there is no prior reason why the
dependency across risks should be the same as the dependency across multiple spells. This is the motivating point of our approach.
This paper contributes as follows. First, it adopts a nested Archimedean copula framework (compare Hofert and M¨achler, 2011; Hofert, 2012) to model risk and multiple spells dependen- cies. Second, it is shown that only the dependency between multiple spells is identifiable and can be estimated. Similar to the single spell models, the competing risks dependence structure is unidentifiable and the presence of multiple spells does not change this. Third, we present a multiple-stage procedure to simulate data and conduct a simulation study to investigate the performance of the proposed model and to assess how misspecification of various model com- ponents lead to estimation bias in quantities of interest. It is shown that estimated risk-specific latent marginals and the partial effects of covariates on these marginals critically depend on the assumed competing risks dependence structure. The MPHM that assumes identical depen- dence structures between risks and spells produces only a snapshot and ignores the variety of possible result patterns.
The paper is structured as follows. Section 2 introduces the model. Section 3 presents the results of a simulation study and Section 4 contains an application to seasonal unemployment duration data. The last section provides additional remarks and ideas for extensions.
2 The Model
We restrict the model to two competing risks r = 1,2. An extension to more than two risks is straightforward by applying risk pooling (Lo and Wilke, 2014). For the k’th spell with k = 1, ..., K there is one pair of latent durations for risks 1 and 2. Let Trk be the latent duration for risk r in the k-th spell. The number of spells is random with Pr(K ≥2)>0 and let k∗ be the realised number of spells. Given that K = k∗ for some k∗ ≥ 1 we observe k∗ covariates x1, ...,xk∗ in IRp. Regressorsxk are allowed to vary across spells but are assumed to not vary within a spell. Let x = (x01, ...,x0K)0. Given that K = k∗ and given the values of x, the conditional distribution of the tuple (T1k, T2k) only depends on xk and not on k∗ and not on (xl : l 6=k). Conditionally on K =k∗ and x, Trk has latent marginal conditional survival function Sr(trk;xk) = Pr(Trk > trk|K =k∗,x) for r= 1,2 and k ≤k∗. Let Tk = min{T1k, T2k} be the minimum of the two latent durations of the k’th spell. The overall survival function
S(t;xk) = Pr(Tk > t|K = k∗,x) = Pr(T1k > t, T2k > t|K = k∗,x) is the conditional survival function of Tk.
Example: Unemployment Duration. Let r = 1 correspond to an exit to employment
and r = 2 to labour force withdrawal. T1k is the latent duration of unemployment until exit to employment in the k’th unemployment spell. The latent duration is not necessarily the observed duration because if unemployment was terminated earlier by an exit to risk 2, T1k is not observed but T2k. Thus, there is a missing data problem as only the shortest latent duration is observed for each spell, which is Tk (observed k’th duration). It is for this reason that Sr(trk;xk) for r = 1,2 is not identifiable unless additional restrictions are imposed. In contrast, the overall survival or survival of the minimum, S(t;xk)is identified for all k ≤k∗.
Let HM(t11, t21, ..., t1k∗, t2k∗;k∗,x) = Pr(Trk > trk for r = 1,2 and k = 1, ..., k∗|K = k∗,x) be the joint conditional survival function of all latent durations Trk and let HD(t1k, t2k;xk) = Pr(Trk > trk for r = 1,2|K = k∗,xk) be the joint survival function of two competing latent durations for the k-th spell and for k ≤ k∗. Note that under our assumptions HD(t1k, t2k;xk) depends neither onK nor on (xl :l6=k). ForHM and HD we assume that they have a nested copula structure:
HM(t11, t21, ..., t1k∗, t2k∗;k∗,x)
=CM[HD(t11, t21;x1), . . . , HD(t1k∗, t2k∗;xk∗)], HD(t1k, t2k;xk) = CD[S1(t1k;xk), S2(t2k;xk)],
where k= 1, ..., k∗ and whereCM andCD are copula functions. We refer to CM as the mother copula and CD as the daughter copula.
Example: Continued. CM characterises any dependencies between theTks for thek∗ spells.
As motivated in the introduction the dependence can be due to omitted variables or frailty that play a role for Tk after conditioning onxk. It should be expected that the dependence is positive. CD characterises the dependence between T1k and T2k given k and xk for all k ≤k∗. This dependence can also be the result of frailty. As explained in the introduction, the effects of frailty can be different for each risk and therefore the direction of the dependence is unknown.
Therefore, it is important to accommodate different dependence structures in CD and CM.
We make the assumption that both copula functions are Archimedean, i.e. Cl(u1, u2) with l ∈ {M, D}fulfills:
Cl(u1, ..., uk) =ψl−1[ψl(u1) +...+ψl(uk)]
wherek =Kork = 2, respectively, and whereψl(u) : [0,1]→IR+is a known strictly decreasing function. ψl is twice differentiable with ψl(1) = 0 and its inverse ψl−1 is completely monotonic on [0,∞].
Under these conditionsHM and HD are survival functions, see Nelsen (2006). Thus, we get the following nested Archimedean copula structure (Hofert, 2012):
HM(t11, t21, ..., t1k∗, t2k∗;k∗,x)
=ψ−1M[ψM(HD(t11, t21;x1)) +. . .+ψM(HD(t1k∗, t2k∗;xk∗))], (1) HD(t1k, t2k;xk) = ψD−1[ψD(S1(t1k;xk)) +ψD(S2(t2k;xk))],
where k = 1, ..., K and where ψM and ψD are the generator functions of the mother and daughter copula, respectively. Since the functions ψM and cψM for a constant c > 0 define the same copula function, we add the norming constraint ψM0 (1) = ψD0 (1) = −1. Figure 1 illustrates the nested copula structure for the 2 risks model with 3 multiple spells.
Figure 1: Nested Archimedean copula structure of a competing 2 risks model with three re- peated occurrences.
CM(·)
CD(·) u11 u12
CD(·) u21 u22
CD(·) u31 u32
TheK overall survival functions for each spell are linked by the so called mother copulaCM. Within each of the multiple spells, the competing risks are linked by the so called daughter copula such that S(tk;xk) = CD[S1(t1k;xk), S2(t2k;xk)]. This model setup implies that the dependence across multiple spells is separate from the dependence between competing risks when CM 6=CD. It therefore differs from other multiple spells models. See e.g. Honor´e (1993), Horny and Picchio (2010) and Lo (2015) who consider models where the dependency between
risks and spells is the same, which correspond to the case of CM =CD =C in our model.
For identifiability we do not make parametric assumptions on Cl. For estimation and the data analysis we will work with parametric copulas. For more details on copulas and properties of Archimedean copulas see for example Nelsen (2006). Several copulas and related equations that are relevant for the subsequent simulation studies are given in section S.I of the supple- mentary material.
The results of this paper carry over to a model with more than two risks without additional assumptions. This is because an R-dimensional Archimedean copula (with R > 2) can be rewritten as a two competing risks model by pooling any R−1 risks into a single risk having no direct interpretation while keeping the last risk that is of interest (Lo and Wilke, 2014).
The copula graphic estimator bases on observable cumulative incidencesQr(t;xk) = Pr(Trk ≤ t, r;xk) withr = 1,2 and assumed CD to reveal S1 and S2. We follow this approach under the additional restriction that CD is Archimedean because then there is a closed form solution for Sr with r= 1,2 (Rivest and Wells, 2001). Carri`ere (1995) considers a more general but known CD and shows that Sr can be obtained by solving a set of nonlinear differential equations. If CD was unknown, these approaches only set identifySr and the bounds are typically very wide.
In order to identify CD, full or partial knowledge of Sr is required (Carri`ere, 1994)). In this paper we follow Zheng and Klein (1995) and leaveS1 andS2 unspecified. We have the following identification result.
Proposition 1 In the described model where one observes K, x, 1T11<T21, ...,1T1K<T2K and T1, ..., TK and where, in particular, (1) holds, the generator function ψM of the mother copula CM is identifiable.
The proof of Proposition 1 is given in Appendix I. The proof only requires CM to be Archimedean, while no restriction on CD is being made. This means ψM is identifiable in a more general nested copula structure with non Archimedean daughter copula.
Put
∆k=
1 if T1k < T2k
2 if T1k ≥T2k and ∆ = (∆1, ...,∆K)T and T = (T1, ..., TK)T.
Proposition 2 The conditional distribution of (T,∆), givenK and x, coincides with the con- ditional distribution of a data generating process, where for k = 1, ..., K the random variables T1k and T2k are conditionally independent, given K and x.
The proof of Proposition 2 is also given in Appendix I. It follows from Proposition 2 that the nonidentifiability of ψD in the single spell model (Tsiatis, 1975) carries over to the multiple spells setting. Proposition 2 also shows that one cannot achieve identification of ψD by making parametric assumptions as long as the parametric model also includes the case thatT1k andT2k are conditionally independent. Thus, the data do not contain information on the dependency betweenT1kand T2k. As already discussed above, our model can be fitted for different assumed dependence structures betweenT1kandT2k. This gives insights how parameters of interest vary with different specifications of dependence.
Although Proposition 1 uses the assumption of a nested Archimedean copula, the iden- tification result also holds for related but more general copula structures. For example, in a hierarchical Archimedean copula, the daughter copula is allowed to be spell-specific and is different from the mother copula. We have
HM(t11, t21, ..., t1k∗, t2k∗;k∗,x)
=CM{CD1[S1(t1k;xk), S2(t2k;xk)], . . . , CDk∗[S1(t1k∗;x∗k), S2(t2k∗;x∗k)]},
where CD1, . . . , CDk∗ are different. The mother copula CM can still be identified from the dependency of observable spells. In contrast, the daughter copulasCDk must be known as usual for the Copula Graphic estimator. This extension allows for more flexibility in the dependence structure, but still requires more knowledge about the daughter copula. This is inevitable due to the fundamental non-identifiability of competing risks model.
Estimation We assume that the copulas are known up to their unknown parameters θl for l ∈ {M, D}. θl is of dimension 1 or 2 for common Archimedean copulas. For one parameter copulas there is a 1-1 relationship between θ and Kendall-τ ∈ [−1,1], which is a measure of pairwise correlation between two random variables that are linked by a copula. While we estimate the parameters of the copula, we normally report Kendall-τ in the simulations and application for reasons of convenience.
Suppose there is a random sample of units i = 1, ..., m. The random number of multiple spells for unit i is ni, i.e. the number of observable spells may differ across i. In typical applications ni is a single-digit or two-digit count variable. In order to be able to estimate θM it is assumed that ni >1 for at least some i. Let l =Pm
i=1ni be the total number of observed spells in the sample. The observed regressors are xik, the stacked matrix of regressors is X which isl×p. This means that the values of regressors are allowed to vary with multiple spells.
Tik is the observable duration for the k’th spell of unit i. δik is the risk indicator for the k’th-spell of thei’th unit and defined byδik = 1 if argminr={1,2}{T1ik, T2ik}= 1 and 0 otherwise for alliandk. The latent durations{T1ik, T2ik}are not observed but onlyTik =min{T1ik, T2ik}.
Another reason for not observing{T1ik, T2ik}is independent censoring, for example, at the end of the observation period. This can be easily incorporated in the usual way but is omitted here for notational convenience.
As usual for the copula graphic estimator, models for Qr and not for Sr are estimated.
Several parametric models for Qr and related equations that are relevant for the subsequent simulation studies and application are given in section S.II of the supplementary material. Qr are known functionals up to unknown nuisance parameters γr ∈ Γ ⊂ IRq and parameters on the regressors βr∈IRp. γ and β are the stacked parameter vectors for all risks.
Let Λr(t) =Rt
0 λr(u)du and λr be the cause-specific hazard function for r = 1,2. Λ(t;x) = Rt
0λ(u;x)duis the integrated hazard function such thatλ=λ1+λ2. The pseudo cause-specific survival function (Jeong and Fine, 2006) is
Wr(t;x) = exp(−Λr(t;x)) (2) such that S(t;x) =W1(t;x)W2(t;x) . The cumulative incidence function forr = 1,2 is
Qr(t;x) = Z t
0
λr(u;x)S(u)du. (3) Qr is related to the copula by
Q0r(tk;x) = −∂CD(S1(t1k;x), S2(t2k;x);θD)
∂trk
t
1k=t2k=tk
. (4)
The joint density function of (Tk, δk) is given by the K-order derivative of (1), f((t1, δ1), ...,(tK, δK))
= ∂KCM(S(t1), ..., S(tK))
∂t1...∂tK
= ∂KCM(CD(S1(t11), S2(t21))), ..., CD(S1(t1K), S2(t2K))
∂CD(S1(t11), S2(t21))...∂CD(S1(t1K), S2(t2K))
CD(S1(t1k),S2(t2k))=S(tk),∀k=1,...,K
× Y
k|δk=1
∂CD(S1(t1k), S2(t2k))
∂t1k
t1k=t2k=tk
Y
k|δk=0
∂CD(S1(t1k), S2(t2k))
∂t2k
t1k=t2k=tk
= ∂KCM(CD(S1(t11), S2(t21))), ..., CD(S1(t1K), S2(t2K))
∂CD(S1(t11), S2(t21))...∂CD(S1(t1K), S2(t2K))
CD(S1(t1k),S2(t2k))=S(t),∀k=1,...,K
×
K
Y
k=1
Λ01(t1k)δkΛ02(t2k)(1−δk)S(tk). (5)
Given the density in (5), the sample likelihood function is
L(θM,γ,β;X)
=
m
Y
i=1
Li(θM,γ,β;X)
=
m
Y
i=1
∂niCM(CD(S1(t1i1;xi1), S2(t2i1;xi1),θD), ...
∂CD(S1(t1i1;xi1), S2(t2i1;xi1),θD), ...
..., CD(S1(t1ini;xini), S2(t2ini;xini),θD),θM) ...∂CD(S1(t1ini;xin1), S2(t2ini;xini),θD)
CD(S1(t1ik;xik),S2(t2ik;xik),θD)=S(tik;xik)
×
ni
Y
k=1
(Λ01(tik;xik))δik(Λ02(tik;xik))(1−δik)S(tik;xik)
#
. (6)
From (6), thei’th unit’s contribution to the log likelihood is
lnLi(θM,γ,β;X) = lnL1i(θM,γ,β;X) + lnL2i(γ,β;X)
fori= 1, ..., m. The log likelihood can be therefore decomposed into two parts. lnL1 is the part of the model related to CM(·). lnL2 is the part of the model related to Qr(t;x) and does not depend on the copula structure. Thus, the likelihood in (6) can be maximised in one step or with a two step procedure. In the latter caseL2 is first maximised to obtain the usual partial or pooled MLE in the context of dependent data. Given the estimated γ andβ,L1 is maximised.
In particular for semiparametric models for Qr, a 2-step procedure is more convenient but the one step procedure is more efficient. In our simulations and application we estimate in one step.
Let ˆθM, ˆγ and ˆβ be the MLE maximising (6). Provided that it is a parametric estimator, the usual asymptotic properties apply when m → ∞ and P(ni > 1) > 0. The MLE is asymptotically normal and consistent with the usual variance matrix. It is remarked that we have fully specified the conditional density (full MLE). The likelihood in (6) therefore reflects the spell dependence and thus no further corrections or adaptations arising from the dependence of spells are required. This is in contrast to the case when one ignores any spell dependencies by using a partial likelihood approach. In the case of independent multiple spells, the copula structure reduces toCM(S(t1), . . . , S(tK)) =CD(S1(t11), S2(t21))∗...∗CD(S1(t1K), S2(t2K)) with lnL1 = 0. This partial likelihood corresponds to a pooled model for the cumulative incidence that treats each spell in the sample as a single spell. Ignoring any spell dependencies in an application therefore invalidates the estimated variances of the coefficients. This is because the scores of the partial likelihood are correlated across spells and the information matrix equality no longer holds. A practical alternative to specifying the full density, as in (6), is to maximise lnL2 only and to compute cluster robust standard errors and inference statistics, where the clusters are formed by the multiple spells from the same unit.
Similar to panel data with the number of periods going to infinity, the MLE of (6) could be also obtained for m fixed and nonrandomni with ni → ∞ for at least onei. In the case of m = 1, the panel analogue would be a time series. However, such a scenario is impractical in applications as it requires very high order derivatives of the copula generator which are difficult to obtain (compare section S.I of the supplementary material). It should be also irrelevant for most applications with survival data as ni is typically a small integer and much smaller than m. In applications with ni → ∞ one may also prefer time series models where dependence between (T1k, T2k) and (T1k0, T2k0) decreases for increasingk0−k > 0.
3 Simulation study
We conduct a comprehensive simulation study to investigate the features of various variants of the model of Section 2. Existing simulation designs for multiple spells models typically only
consider the cause specific (e.g. Sankaran and Anisha, 2011) or sub- distributions (e.g. Huang et al., 2017) with dependencies that come from parametric frailty distributions. Simulating data for our suggested model is complicated as it requires known cumulative incidences and marginals under an assumed nested copula structure. The underlying nested Archimedean copula model is difficult to simulate as the dependence structure mixes two different copulas. We simulate data using a multiple step procedure that is outlined in detail in section S.III of the supplementary material. We base on the procedure of Lo (2015) for a common spell and risk copula. The mixture of the copulas is simulated using the framework for sampling nested Archimedean copulas outlined by Hofert and M¨achler (2011) and Hofert (2012). For simplification we only simulate models in which the copula generator of the mother and daughter copula coming from the same family (but possibly with different parameters) as otherwise the copula mixture becomes intractable. For the Gumbel copula we simulate data with different combinations of τM and τD. For tractability we restrict to the case τM = τD for the other copulas. Since we mainly focus on investigating the dependence structure and the implied latent marginal survivals, we focus on models without regressors. We limit our analysis to models with up to three multiple spells, i.e. ni ≤3 for alli, due to the complexity of the higher order derivatives of the inverse of the copula generator (compare section S.I of the supplementary material).
In our models, unless otherwise stated, we draw data for m = 200, where 50% have n1 = 1, 25% have ni = 2 and the rest hasni = 3. ni is random for alli. We therefore have l= 350 and m/l = 57%. In all simulation designs we sample 500 times.
We simulate data for several popular copulas and models for the cumulative incidences.
We investigate finite sample performance and the role of the ratio of number of multiple spells and the number of units. Beside looking at correctly specified models, we study the role of misspecification of the dependence structure and the model for the cumulative incidence. In particular, we use the Frank copula, the Clayton copula and the Gumbel copula. The models forQrcomprise the log-normal accelerated failure time model, the log-logistic proportional odds model and the odd-rate transformation model with Gompertz cause specific baseline hazard.
Relevant details of the copulas are given in section S.I of the supplementary material and those of the models for Qr are given in section S.II of the supplementary material.
We assess the performance of the estimates forτM andSrby computing their bias and mean squared error (MSE). We focus on the properties of the implied Sr and not the estimation of
the parameters of the model forQr as only the former are affected by the (in-)correct modelling of the dependence structures. While the bias and the MSE are obvious to compute forτM, the measures forSr vary int. To simplify the reporting we display their averages overt which have a global character. In particular, the average bias (AB) is
AB(Sr) = 1 J
J
X
j=1
EˆSˆr(tj)−Sr(tj)
and the average mean squared error (AMSE) is
AM SE(Sr) = 1 J
J
X
j=1
E[( ˆˆ Sr(tj)−Sr(tj))2]
for risks r = 1,2, where tj = t1, . . . , tJ are equidistant grid points on the support of T with J=100. We estimate the expected values by taking the average over the 500 simulations. As the Sr partly vary across models due to different DGPs, the Bias and MSE should be related to their magnitude for comparability. We therefore also consider a relative AMSE (RAMSE) which has a percentage interpretation. It is obtained by dividing the AMSE(Sr) through S¯r = PJ
j=1Sr(tj). It is zero if the AMSE was zero and it is one if the AMSE had the same value as ¯Sr.
Table 1 shows the results for 9 different models. Some of them are correctly specified (Models 1 and 2), some have a misspecified τD (Models 3-5), some have a misspecified copula generator (Models 6,7), a misspecified model for Qr (Model 8) and a misspecified ψ and Qr (Model 9). Models 3-9 therefore assess how important misspecification of the various model components is. Model 3 is of particular interest as it makes the same assumption as commonly used duration models that τD is unknown but the same as τM. The results suggest that the model is well estimated for various degrees of dependencies, various copulas and various models for Qr when they are correctly specified (Models 1-2 and compare the additional results in section S.IV of the supplementary material). Misspecification of τD leads to by far the largest bias and AMSE, while misspecification of ψD,ψM orQr only leads to smaller bias and smaller increases in AMSE. Imposing the restriction τM =τD (Model 3) can therefore lead to a sizable bias if the two differed in the true model. If the model is estimated under this restriction without knowing the true τD, it produces only a snapshot. The uncertainty of not knowing it
is not reflected.
There is a discussion in the applied econometrics literature on the matter of omitting single spells in a multiple spells model and how the presence of multiple spells contributes to model identification. For this reason the results for the same models but based on less data are shown in Table 2. These results are obtained by only using the first two spells of each unit. This means that we omit those units with one spell only and omit the 3rd spell of a unit, if any. The results suggest that, when compared to Table 1, the AMSE becomes larger but, in most of the time, only slightly, which thus confirms that omitting relevant information decreases precision of statistical results. But it does not bias the results as long as the omission is random (i.e.
not related with T1ik and T2ik for k = 1, ..., K) as in our simulation design.
Table 3 analyses two questions. First, to what extent it matters for the precision of results when the share of single spells (units with ni = 1) in the sample increases, keeping the total number of spells (l) in the sample constant. This is of interest as single spell units do not contribute to the estimation of τM and a small share of multiple spells may therefore affect the precision of all estimates. Second, we relate the estimates from the multiple spells model to estimates that are obtained under the assumption that all spells are independent. In this case τM = 0 is assumed and not estimated. We simulate two models with m/l = 34% and 94%.
The results show that the MSE for ˆτM increases strongly (six times larger) when the number of multiple spells in the sample drops. In contrast, the MSE for ˆS1 and ˆS2 increases by factor 3 when the share of single spells drops. The table also contains results when a single spell model is estimated that assumes all spells are independent. For the model that assumes independent spells the MSE for ˆS1 and ˆS2 is only marginally higher than for the multiple spells model.
This suggests that the main benefit of employing a multiple spells model is that it produces an estimate for the degree of dependence between spells.
Table1:BiasandMeanSquaredErrorrelatedmeasuresforvariousmodelswithm=200andl=350.Assumedmodelfeatures correctunlessotherwisestated. Model123 DGP:Gumbel,LNAFT,γ1=(2,0.5),γ2=(1.5,2),τM=0.3 DGPτD=0.3τD=0.7τD=0.7 AssumedcorrectcorrectτM=τD (A)B(A)MSERAMSE(A)B(A)MSERAMSE(A)B(A)MSERAMSE ˆτM-0.00260.0021-0.00010.0023-0.00010.0023 ˆS10.00110.00010.19%-0.00020.00010.14%0.03410.00244.82% ˆS20.00020.00030.30%-1.62e-050.00010.17%0.05720.006912.48% Model456 DGP:Gumbel,LNAFT,γ1=(2,0.5),γ2=(1.5,2),τM=0.3 DGPτD=0.7τD=0.7τD=0.3 AssumedτD=0.3τD=0.0Frank (A)B(A)MSERAMSE(A)B(A)MSERAMSE(A)B(A)MSERAMSE ˆτM-0.00010.0023-0.00010.00230.00050.0024 ˆS10.03410.00234.76%0.08380.011423.18%0.00060.00020.25% ˆS20.005690.006712.19%0.13780.029754.01%0.00130.00040.40% Model789 DGP:Gumbel,LNAFT,γ1=(2,0.5),γ2=(1.5,2),τM=τD=0.3 AssumedClaytonORGOMFrank,LLPOM (A)B(A)MSERAMSE(A)B(A)MSERAMSE(A)B(A)MSERAMSE ˆτM-0.09930.01190.01660.0023-0.00680.0024 ˆS1-0.02440.00091.22%0.00430.00050.60%0.00380.00030.42% ˆS2-0.03530.00212.18%0.00130.00080.64%-0.00540.00050.52%
Table2:BiasandMeanSquaredErrorrelatedmeasuresforthemodelsinTable1onlyfor1stand2ndspellforunitswithni>1, i.e.m=100andl=250.Assumedmodelfeaturescorrectunlessotherwisestated. Model123 DGP:Gumbel,LNAFT,γ1=(2,0.5),γ2=(1.5,2),τM=0.3 DGPτD=0.3τD=0.7τD=0.7 AssumedcorrectcorrectτM=τD (A)B(A)MSERAMSE(A)B(A)MSERAMSE(A)B(A)MSERAMSE ˆτM-0.00770.0036-0.00140.0037-0.00140.0037 ˆS10.00070.00030.34%-0.00060.00010.24%0.03530.00254.55% ˆS2-0.00010.00060.52%-0.00030.00020.28%0.06250.007912.78% Model456 DGP:Gumbel,LNAFT,γ1=(2,0.5),γ2=(1.5,2),τM=0.3 DGPτD=0.7τD=0.7τD=0.3 AssumedτD=0.3τD=0.0Frank (A)B(A)MSERAMSE(A)B(A)MSERAMSE(A)B(A)MSERAMSE ˆτM-0.00140.0037-0.00140.0037-0.00380.0040 ˆS10.03510.00244.45%0.08550.011521.11%-0.00030.00030.38% ˆS20.06180.007612.32%0.14730.032853.46%-3.95e-050.00070.58% Model789 DGP:Gumbel,LNAFT,γ1=(2,0.5),γ2=(1.5,2),τM=τD=0.3 AssumedClaytonORGOMFrank,LLPOM (A)B(A)MSERAMSE(A)B(A)MSERAMSE(A)B(A)MSERAMSE ˆτM-0.10960.01540.01260.0034-0.00930.0037 ˆS1-0.02570.00111.25%0.00380.00060.68%0.00390.00030.44% ˆS2-0.03920.00262.29%0.00310.00090.83%-0.00560.00050.56%
Table 3: Bias and Mean Squared Error related measures when the number of multiple spells in the sample decreases while the total number of spells remains constant at l=350. Assumed model features correct.
DGP: Gumbel, LNAFT, γ1 = (2,0.5), γ2 = (1.5,2), τM =τD = 0.7
m 119 330
#ni 3/1/115 317/6/7
m/l 34% 94%
(A)B (A)MSE RAMSE (A)B (A)MSE RAMSE
ˆ
τM -0.7642 0.5456 3.9881 3.2370
Sˆ1 -0.2347 0.1525 0.28% 0.6248 0.0651 0.12%
Sˆ2 -0.0887 0.1934 0.32% 0.7206 0.0830 0.13%
Assumed single spell model
Sˆ1 -0.1845 0.1582 0.29% 0.6230 0.0698 0.13%
Sˆ2 -0.1181 0.1979 0.33% 0.7383 0.0881 0.14%
Note: Values for (A)B and (A)MSE are in ×10−3.
4 Analysis of multiple spells unemployment
In this section we put the model to real data. The empirical analysis of unemployment duration is of great societal importance which is documented by a large amount of academic publications.
While most analysis relies on single spell data, a number of contributions have used multiple spells data (for example Rød and Westlie, 2012, Kalwij 2010, Carrasco and Garcia-Perez, 2015).
It is therefore interesting to see what insights the suggested models of Section 2 produce.
In our application we consider unemployment duration of seasonally employed workers in Germany. Seasonal employment is common in the German labour market and a number of policy changes have been undertaken to mitigate winter unemployment due to temporary lay- offs (Arntz and Wilke, 2012). We extend previous analysis by estimating a competing risks model.
For the empirical analysis we use the scientific use file version of the sample of the inte- grated labour market biographies (SIAB) 1975-2010 of the Institute for Employment Research (IAB), Germany. These administrative data are a 2% random sample of the workforce in Ger- many that contributed to the social insurance in the period 1975-2010. The SIAB contains daily information about periods of dependent employment and claim periods for unemploy- ment compensation along with basic information about the individual (such as gender, wage, age and employment history) and the employing firm (such as business sector and location).
For more information on the SIAB see Dorner et al. (2010). From these data we extract all unemployment benefit claim periods. We define unemployment benefit duration as receipt of unemployment benefits (Arbeitslosengeld) from the German Federal Employment Agency. We only use observations that were generated by aged 16-65 in the most northern federal states of Hamburg or Schleswig-Holstein. We also only consider unemployment benefit claim periods if the previous employment period had been in a business sector prone to seasonal employment patterns (construction and agriculture) and if the previous employment period ended during the winter season. After some data preparation steps we are left with a sample of 2,182 individ- uals which start in total 3,786 1st, 2nd and 3rd spells during the winter seasons in the period 1981-2009. Thus, we omit higher order spells from the analysis as the likelihood becomes oth- erwise untractable. There are 2 right censored spells that we recode to unknown for simplicity.
Similar to the simulation study we work with 3 different estimation samples. One that only comprises the individuals with 2 spells, one that omits the third spells and one that uses all 1st-3rd spells. We define the following 5 exit states (competing risks) in the data:
1. New seasonal employer (in construction or agriculture) 2. Recall (to the pre unemployment employer)
3. Other business sector
4. Other benefits (unemployment assistance, training measures) 5. other/unoberved (no data available)
The model is estimated 4 times (for risks 1-4) with 2 risks as the so called risk pooling approach (compare Lo and Wilke, 2014) can be applied without additional restrictions. We do not consider explicit estimates for the 5th risk as it is a pooled remainder state and does not have a clear interpretation. We estimate the model with assumed Gumbel copula and Log-normal accelerated failure time model forQr. The resulting parameter estimates are reported in Table 4.
The results suggest that there is a small but significant positive rank correlation between spells (ˆτM ≈ 0.04−0.08). This can be interpreted as that the length of spells is only weakly correlated across multiple occurrences. This may appear a bit surprising but can be explained
Table 4: Estimated coefficients of seasonal unemployment duration model.
only 2 spells 1st & 2nd spells 1st-3rd spells indep. spells
m 1,045 2,182 2,182
#ni = 1/2/3 0/1,045/0 2,182/1,045/0 2,182/1,045/559
l 2,090 3,227 3,786 3,786
coef. S.E. coef. S.E. coef. S.E. coef. S.E.
risk 1 ˆ
γ11 1.301∗∗∗ 0.043 1.481∗∗∗ 0.042 1.479∗∗∗ 0.040 1.461∗∗∗ 0.038 ˆ
γ21 6.374∗∗∗ 0.252 6.349∗∗∗ 0.219 6.213∗∗∗ 0.199 6.220∗∗∗ 0.194 βˆ11 (female) 1.138∗∗∗ 0.269 1.335∗∗∗ 0.210 1.211∗∗∗ 0.188 1.182∗∗∗ 0.186 βˆ21 (German) -0.585∗∗ 0.243 -0.281 0.211 -0.137 0.192 -0.155 0.187 βˆ31 (educ) -0.219∗∗ 0.097 -0.223∗∗ 0.089 -0.239∗∗∗ 0.083 -0.242∗∗∗ 0.081 θˆM 1.053∗∗∗ 0.018 1.072∗∗∗ 0.023 1.082∗∗∗ 0.020
ˆ
τM 0.051 0.067 0.076
lnL -12,819.119 -20,141.848 -23,553.887 -23,571.358
risk 2 ˆ
γ12 1.046∗∗∗ 0.025 1.222∗∗∗ 0.027 1.189∗∗∗ 0.024 1.187∗∗∗ 0.023 ˆ
γ22 4.973∗∗∗ 0.131 5.175∗∗∗ 0.131 5.122∗∗∗ 0.119 5.146∗∗∗ 0.117 βˆ12 (female) -0.342∗∗∗ 0.114 -0.002 0.100 -0.030 0.090 -0.064 0.088
βˆ22 (German) 0.016 0.131 0.158 0.131 0.126 0.120 0.108 0.117
βˆ32 (educ) -0.010 0.063 -0.019 0.063 -0.023 0.056 -0.033 0.055 θˆM 1.052∗∗∗ 0.018 1.070∗∗∗ 0.022 1.079∗∗∗ 0.019
ˆ
τM 0.049 0.065 0.073
lnL -13,122.865 -20,611.272 -24,145.436 -24,162.167
risk 3 ˆ
γ13 1.637∗∗∗ 0.080 1.655∗∗∗ 0.053 1.699∗∗∗ 0.054 1.673∗∗∗ 0.051 ˆ
γ23 6.338∗∗∗ 0.293 6.125∗∗∗ 0.218 6.250∗∗∗ 0.217 6.272∗∗∗ 0.212 βˆ13 (female) 0.508∗ 0.304 -0.153 0.159 -0.060 0.158 -0.114 0.155
βˆ23 (German) 0.256 0.277 0.158 0.211 0.146 0.210 0.111 0.204
βˆ33 (educ) 0.034 0.139 0.128 0.100 0.116 0.098 0.105 0.095
θˆM 1.059∗∗∗ 0.019 1.075∗∗∗ 0.023 1.086∗∗∗ 0.020 ˆ
τM 0.055 0.070 0.079
lnL -12,458.326 -19,958.003 -23,248.979 -23,267.116
risk 4 ˆ
γ14 1.176∗∗∗ 0.055 1.284∗∗∗ 0.045 1.257∗∗∗ 0.041 1.256∗∗∗ 0.041 ˆ
γ24 5.302∗∗∗ 0.186 5.414∗∗∗ 0.161 5.443∗∗∗ 0.150 5.459∗∗∗ 0.149 βˆ14 (female) 0.724∗∗ 0.286 0.088 0.153 0.175 0.146 0.167 0.146 βˆ24 (German) 0.855∗∗∗ 0.181 0.895∗∗∗ 0.157 0.858∗∗∗ 0.147 0.843∗∗∗ 0.145 βˆ34 (educ) 0.239∗∗ 0.107 0.222∗∗ 0.090 0.188∗∗ 0.083 0.183∗∗ 0.082 θˆM 1.046∗∗∗ 0.017 1.062∗∗∗ 0.021 1.071∗∗∗ 0.018
ˆ
τM 0.044 0.058 0.067
lnL -12,344.408 -19,585.954 -22,876.09 -22,891.36
by the construction of our sample of seasonally unemployed and by studying the length of unemployment benefit claim spells. The maximum entitlement length for benefits is legally restricted and depends on the length of employment periods before unemployment and on age.
It is also evident from Table 4 that coefficients change across risks, justifying the use of a competing risks model. The results also confirm that only using individuals with at least 2 spells causes a sample selection problem, i.e. dropping those with only one spell. In particular, the coefficient on female appears to be considerably biased compared to when individuals with only one spell are also included. The reason for the sample selection bias is that if the first spell was long, an individual will less likely produce a second spell in the observation period. For this reason, dropping those with one spell is related to the dependent variable and therefore not random conditional on the regressors. Including the 3rd spells in the analysis makes results a bit sharper but generally they remain similar. Treating all spells as independent slightly underestimates the standard errors. The more positive τM, the more are the standard errors expected to be underestimated when the dependency between multiple spells is ignored.
In what follows we only focus on the estimates obtained from the sample with all 1st-3rd spells as it appears to produce the most precise estimates. We compute implied ˆSr(t;x = (0,0,0), τd) which are shown in Figure 2. It is apparent that the assumed τD play a big role for the resulting estimates, confirming the findings of the simulation section and for single spell models. Thus, there is little evidence that the use of multiple spells helps in mitigating the non-identifiability of the competing risks model.
To illustrate the effect of assuming τD on partial effects of regressors on Sr, we compare Sˆr(t;x1) with ˆSr(t;x0), i.e.
∆ ˆSr(t;τD) = ˆSr(t;x1, τD)−Sˆr(t;x0, τD),
for x1 = (1,0,0) and x0 = (0,0,0). This means we compare the estimated latent marginal conditional survival curves for females with the one of the males holding the other variables at 0. The resulting partial effects are shown in Figure 3. It is evident that the assumed τD plays an important role for the effects as the magnitude and the sign depend on it. Thus, not only the distribution but also partial effects may be misleading when only one value of τD, such as ˆ
τM is considered. In order to summarise the estimated partial effects in one figure we compute
Figure2:EstimatedˆSr(t;x=(0,0,0),τD)astheyvarywiththeassumedτD.
0 .2 .4 .6 .8 1
0200400600800 unemployment duration (in days) S1(t;x=0,τD=.01)S1(t;x=0,τD=.1) S1(t;x=0,τD=.2)S1(t;x=0,τD=.3) S1(t;x=0,τD=.4)S1(t;x=0,τD=.5) S1(t;x=0,τD=.6)S1(t;x=0,τD=.7) S1(t;x=0,τD=.8)S1(t;x=0,τD=.9) S1(t;x=0,τD=τM)
New seasonal employer
0 .2 .4 .6 .8 1
0200400600800 unemployment duration (in days) S2(t;x=0,τD=.01)S2(t;x=0,τD=.1) S2(t;x=0,τD=.2)S2(t;x=0,τD=.3) S2(t;x=0,τD=.4)S2(t;x=0,τD=.5) S2(t;x=0,τD=.6)S2(t;x=0,τD=.7) S2(t;x=0,τD=.8)S2(t;x=0,τD=.9) S2(t;x=0,τD=τM)
Recall
0 .2 .4 .6 .8 1
0200400600800 unemployment duration (in days) S3(t;x=0,τD=.01)S3(t;x=0,τD=.1) S3(t;x=0,τD=.2)S3(t;x=0,τD=.3) S3(t;x=0,τD=.4)S3(t;x=0,τD=.5) S3(t;x=0,τD=.6)S3(t;x=0,τD=.7) S3(t;x=0,τD=.8)S3(t;x=0,τD=.9) S3(t;x=0,τD=τM)
Other Business Sector
0 .2 .4 .6 .8 1
0200400600800 unemployment duration (in days) S4(t;x=0,τD=.01)S4(t;x=0,τD=.1) S4(t;x=0,τD=.2)S4(t;x=0,τD=.3) S4(t;x=0,τD=.4)S4(t;x=0,τD=.5) S4(t;x=0,τD=.6)S4(t;x=0,τD=.7) S4(t;x=0,τD=.8)S4(t;x=0,τD=.9) S4(t;x=0,τD=τM)
Other benefits
