Joint Modeling of Multivariate Survival Data With an Application to Retirement

Abstract

The Cox proportional hazards model has been pervasively used in many social science areas to examine the effects of covariates on timing to an event. The standard Cox model is intended to study univariate survival data where there is a singular event of interest, which can only be experienced once. However, we may additionally wish to explore a number of other complexities that are prevalent in survival data. For example, an individual may experience events of the same type more than once or may experience multiple types of events. This study introduces innovations in recurrent (repeatable) event analysis, jointly modeling several endogenous survival processes. As an example and an application, we simultaneously model two types of recurrent events in the presence of a dependent terminal event. This model not only correctly handles different types of recurrent events but also explicitly estimates the direction and magnitude of relationships between recurrences and survival. This article concludes with an example of the model to examine how the timing of retirement is associated with the risks of multiple spells of employment and childbearing. The theoretical discussions and empirical analyses suggest that the multivariate joint models have much to offer to a wide variety of substantive research areas.

Keywords

frailty model joint modeling recurrence survival analysis

In the past several decades, a growing recognition of survival analysis has been propelled by the availability of longitudinal or retrospective surveys as well as advances in statistical science. The Cox (1972) proportional hazards model, a widely used procedure in survival analysis, was originally developed by Cox studying the occurrence and timing of terminal events. The issues of timing and sequence of life events are highly relevant for social science researchers and these methods, which incorporate covariate effects on the risk of observing an event are used, for example, to model union formation and dissolution (Pessin 2018; Schneider 2011; Schimmele and Wu 2016), fertility behavior (Axinn, Dirgha, and Smith-Greenaway 2017; Balbo and Barban 2014), and job mobility (Blossfeld and Drobnic 2001; Blossfeld, Hamerle, and Mayer 2014), to name a few.

The standard Cox model applies to univariate survival data when there is a singular event of interest and the event can only be experienced once. In addition, event times are assumed to be statistically independent (Ezell, Land, and Cohen 2003; Therneau and Grambsch 2013). Grappled with increasingly dynamic and heterogeneous life trajectories, social science researchers may additionally wish to explore a number of relevant issues found in multivariate survival data. For example, how does one employ survival analysis to study repeatable events (a singular event type happening more than once) and different types of correlated events?

Unlike many biostatistical studies, in which the event under study may occur only once (e.g., the death of a patient), the majority of events in social sciences are repeatable, such as marriages, divorces, and employments. Recurrent (repeatable) data arise when the same type of event can occur to a subject multiple times. The standard Cox model is not suitable because it assumes a subject is not at risk any longer after the subject experienced the first event. In the analysis of recurrent data, all subjects are at risk of new events as long as they are not censored or have experienced a terminal event, like death (Commenges and Jacqmin-Gadda 2015). In addition, when a subject experiences the same type of event multiple times, the timing of these events is likely to be correlated within subject. Failing to account for the repeatability of the event is tantamount to imposing an independence assumption on the occurrence of the events, which often leads to biased estimation (Box-Steffensmeier and Jones 2004).

Besides repeatable events, there is considerable interest in studying the timing and sequencing of interdependent life events. The outcome of one life event can influence the occurrence of another life event. For example, Lillard’s (1993) research shows that the risk of marital conception is associated with a decreased risk of marital disruption. On the one hand, it is well known that children affect the chances that their parents will divorce (Lillard and Waite 1993). Having children together raises the costs of divorce and increases the gains from marriage, leading to greater marital stability among couples with children (Lillard and Waite 1993). On the other hand, it is also plausible that people take into account anticipated changes in marriage duration in their childbearing decision (Lillard 1993). Therefore, the decision about childbearing and the decision to remain in marriage (or end a marriage) may be subject to shared unobserved factors (Steele 2011). Childbearing outcomes may be endogenous with respect to marital disruption. Estimation of the impact of anticipated events on current event transitions is challenging as these anticipation factors are unobserved (Ermisch and Steele 2016). Petersen (1995) proposed a framework for analyzing multiple types of events where the events are interdependent—that is, where the occurrence of one type of event affects the probability that the other type of event happens. This line of inquiry has greatly extended prior research by developing dynamic modeling for interdependent life events.

Fortunately, in recent years, joint modeling of several survival processes has received considerable attention in statistical research because it can be used to address some interesting scientific questions that could not be answered before, such as the impact of anticipated events on current transitions and quantifying the relationships among them. Separate modeling of each survival process may not fully reveal potential mechanisms and can produce misleading results. Appropriate statistical methods are needed to utilize the richness of these data in order to identify potential relationships between different endogenous survival processes.

Given the ubiquitous repeatable events in social science research, as well as long-standing interests in examining interdependence of life events, we introduce a new dynamic approach to model survival data, which extends the Cox model by incorporating repeated and/or simultaneous events. Most of existing studies using this type of model were conducted in biomedical fields. As a result, these methodologically advanced techniques are relatively unfamiliar to social science researchers. Considering the complexity of dynamic life events, we would like to extend social scientists’ statistical tool kit beyond the standard Cox model by taking advantage of recent developments in multivariate survival models.

Specifically, in this article, we introduce a multivariate joint frailty model for two types of recurrent events in the presence of a dependent terminal event, with right censored survival data. Our application is based on retrospective data from the 2011 Canadian General Social Survey, Cycle 25 (GSS-25) conducted by Statistics Canada. We treat employment and fertility as two repeatable events and discuss how to appropriately handle them using shared frailty models. We aim to assess how the hazards of recurrent events (employment and fertility) would impact the hazard of terminal event (retirement) and at the same time explicitly estimate the strength of associations among these three survival processes. Employment, fertility, and retirement transitions are all specified in separate equations but are estimated in a joint maximum penalized likelihood procedure (see Mazroui et al. 2013). This allows us to analyze the dependencies of the transitions explicitly, controlling for the potential endogeneity of each transition with respect to all the others.

We begin with a review of recurrent events analysis and discuss the use of a shared frailty model to analyze recurrent events. Our attention then turns to an introduction to multivariate joint frailty models and an estimation method in The Multivariate Joint Frailty Model section. Then, we apply multivariate joint frailty models to an analysis of a Canadian national survey with observations of two recurrent events (employment and fertility) and terminal event (retirement) in The Influence of Employment and Fertility on Retirement Timing section. The final section provides conclusions and discussions.

Literature Review

Cox Model

The Cox model analyzes effects of covariates on the hazard rate. Let $h (t; Z)$ be the hazard rate at time t for an individual with risk vector $Z$ . A Cox proportional hazards model is defined as follows:

h (t; Z) = h_{0} (t) exp (β^{T} Z),

where $h_{0} (t)$ is the baseline hazard rate or reference value, $β$ is a parametric vector, and $exp (β^{T} Z)$ is the relative risk, a proportionate increase or reduction in risk. The model is called semiparametric because it does not make specific assumptions about the baseline hazard function.

Parameter estimates in the Cox model are obtained by maximizing the partial likelihood as opposed to the full likelihood. The partial likelihood allows estimation of covariates without making any assumption on baseline hazard. This is a key reason for the popularity of the Cox model. The partial likelihood function is derived by taking the product of the conditional probability of a failure at each time, given the number of subjects that are in the risk set at that time (Cox 1972; Klein and Moeschberger 2005).

Recurrent (Repeatable) Events

In many studies, individuals may undergo the same type of event several times during the follow-up period (see, e.g., Cook and Lawless 2007). Common examples of recurrent events include heart failure hospitalization, asthma attack, marriage, divorce, and unemployment.

Analyzing repeatable events of the same type tends to be more complicated and also raises a number of difficult statistical questions. First, it is necessary to take into account the timing and order of events for the same subject. If one assumes that the first event is no different from the following events, then one may miss important and useful information regarding the timing and sequence of the repeated event (Box-Steffensmeier and Jones 2004). Second, the dependency between these recurrent events should be taken into consideration. As a subject experiences the same type of event more than once, the events from the same subject are potentially correlated (Amorim and Cai 2015; Wienke 2010). The occurrence of one event may change the probability of subsequent events of the same type. This means that the follow-up recurrent events are related to the occurrence of previous events.

The standard Cox model is not suitable for analyzing recurrent events because all events occurring after the first are neglected in the standard Cox model. It would be an inefficient use of data if we only make use of time to first event, ignoring subsequent events (Amorim and Cai 2015). In the analysis of recurrent data, all subjects are at risk of new events as long as they are not censored or have experienced a terminal event, while in the conventional survival analysis, individuals are not at risk after a first event (Commenges and Jacqmin-Gadda 2015).

There has been a recent surge of interest in modeling recurrent events in biomedical research. Marginal models and frailty models have been proposed by statisticians to account for recurrent events in survival analysis. Marginal models are appropriate when the substantive focus is on the effects of covariates rather than the precise nature of the dependence structure as the association between events is considered as a nuisance parameter (Ezell et al. 2003; L. Liu, Wolfe, and Huang 2004). However, in many real-life circumstances, quantifying the dependence structure is the primary focus of research. Frailty models were developed to model the dependence structure of repeated events and can also be easily incorporated into a joint modeling framework of multiple survival process, which will be discussed later in this article.

Frailty Models

Frailty

Demographers James Vaupel, Manton, and Stallard (1979) introduced the concept of frailty and applied it to the study of population mortality. They illustrated that the population hazard does not truly reflect the hazard of individuals from that population. They observed that at the oldest ages, mortality rates show a slower increase even though the hazards for individuals continue to increase. They explained this observed mortality rate decreases as a consequence of the failure-prone—more “frail” —individuals who die at younger ages leaving a subgroup of robust individuals at the older ages. In other nonhealth contexts, we can think of frailty as some observations having a higher risk of experiencing the event than others (i.e., they are more “frail”), but the source of that difference is unobserved. This line of research stimulated a growing body of research that gave priority to the concept of frailty with far-reaching implications (Myrskylä and Fenelon 2012; Palloni and Beltrán-Sánchez 2017; Vaupel et al. 1998; Wrigley-Field 2014).

Univariate frailty models

In the context of survival analysis, frailty models are extensions of standard Cox model. To be more specific, frailty models are multiplicative hazard models consisting of three components: a frailty variable U (random effect), a baseline hazard function $h_{0} (t)$ (parametric or nonparametric), and a term modeling the influence of observed covariates $exp (β^{T} Z)$ (measured covariates; Wienke 2010). The univariate frailty model is defined as follows:

h (t; Z) = h_{0} (t) exp (β^{T} Z + u) = U h_{0} (t) exp (β^{T} Z) .

The frailty $U = exp (u)$ is a nonnegative random variable varying over the population. Frailty distributions are standardized to $E U = 1$ (Wienke 2010). The variance $σ^{2} = Var (U)$ is interpreted as a measure of heterogeneity across the population. When $σ^{2}$ is small, the values of U are closely located around one. When $σ^{2}$ is large, then values of U are more dispersed. Apart from the frailty variable U, all individuals are assumed to follow the same mortality pattern (Wienke 2010).

Shared frailty models

The role of frailty models in recurrent event analysis has received growing attention, where the correlation among event times is a focus of inquiry. The shared frailty model provides an efficient way to model this correlation by introducing a nonnegative frailty variable, U, in the Cox model. The introduced frailty is considered to be shared among the events within the same subject to induce the dependence among them. The dependence is measured through unobservables that are correlated. All survival times that are related to each other have the same level of frailty. Conditional on the frailty, the event times within the same subject are assumed to be independent.

A shared frailty model in survival analysis is defined as follows. Suppose there are n independent individuals and that individual i has n_i observations and associates with the unobserved frailty u_i $(1 \leq i \leq n)$ . The vector $Z_{i j}$ $(1 \leq i \leq n,1 \leq j \leq n_{i})$ contains the covariate information of the event time $T_{i j}$ of the $j th$ observation for the $i th$ individual. Conditional on the frailty term u_i , the survival times in individual i. $(1 \leq i \leq n)$ are assumed to be independent and their hazard functions to be of the form:

r (t | u_{i}) = r_{0} (t) exp (β Z_{i j} + l n (u_{i})) = u_{i} r_{0} (t) exp (β^{T} Z_{i j}) .

The frailties u_i $(i = 1, ..., n)$ are assumed to be independently and identically distributed random variables following some distribution. Hougaard (1995) discussed the choice of frailty distribution for the shared frailty model and noted that the gamma distribution has typically been used to fit the frailty random effect mainly due to mathematical reasons. That is, if one chooses a gamma distribution for the frailty, parameter estimates are easily obtained through likelihood estimation, which is also readily available in R packages version 3.5.3.

Joint modeling framework

Different life events are potentially correlated via dynamic processes. People often take into account anticipated changes in one life event in making decisions in other life events. For example, union transition and childbearing within that union are two related dynamic processes. The decision to end a cohabitation or to move from cohabitation to marriage is likely to be jointly determined with the decision to have a child with that partner (Lillard and Waite 1993). Women might make greater investments in their relationship if they believe they will get married (Lillard and Waite 1993). Estimation of the impact of the anticipated events on current transitions is always challenging because these survival processes are endogenous. Observed and unobserved factors play roles in different event transitions. If decisions about these life events are jointly determined, then we might expect correlated unobserved factors of the models for each process. Joint modeling will serve the purpose as we can explicitly quantify the direction and magnitude of dependencies among these endogenous survival processes.

Joint frailty models

Joint analysis of recurrent event and survival time has received some research attention in recent years as recurrent events are often subject to dependent termination, which has a nonneglectable impact on the occurrences of the recurrent events. For example, the risk of recurrent strokes is positively associated with the risk of death (Lin et al. 2017). To resolve this issue, L. Liu and colleagues (2004) proposed a joint shared frailty model and an estimation method based on the Monte Carlo Expectation-Maximization (EM) algorithm. Along this line of inquiry, Rondeau et al. (2007) proposed jointly modeling the recurrent events and a terminal event using penalized likelihood estimation. They examined episodic relapses of follicular lymphoma (FL) and death, and they found a positive association between these two processes. In other words, those subjects experiencing FL relapses at the highest rate are typically observed for shorter periods of observation due to mortality. By this rationale, they considered the FL relapses and the terminal event process jointly in a joint frailty model setting. It was shown in this article that ignoring the dependence between the recurrence events and the terminal event could lead to unreliable estimates, with regression coefficients falsely nonsignificant or with an underestimation of the recurrence risks.

Following Rondeau et al. (2007), individual correlation between recurrent events and a terminal event is achieved by a shared frailty term. The model can be specified by the hazard functions:

{\begin{array}{l} r_{i} (t | u_{i}) = r_{0} (t) exp (β_{1}^{T} Z_{i 1} + l n (u_{i})) = u_{i} r_{0} (t) exp (β_{1}^{T} Z_{i 1}) & (recurrent event) \\ λ_{i} (t | u_{i}) = λ_{0} (t) exp (β_{2}^{T} Z_{i 2} + γ l n (u_{i})) = u_{i}^{γ} λ_{0} (t) exp (β_{2}^{T} Z_{i 2}) & (terminal event) \end{array},

where $r_{i} (\cdot)$ denotes the hazard of the recurrent events and $λ_{i} (\cdot)$ denotes the hazard of the terminal event. A shared frailty term u_i links the two survival processes together. The random effects u_i are assumed independent and follow a gamma distribution with unit mean and variance $θ$ .

The frailty term acts differently for the two hazard rates (u_i for the recurrent rate and $u_{i}^{γ}$ for the terminal event rate). When $γ = 0$ , the terminal event rate is independent of the recurrent event rate. When $γ = 1$ , the effect of the frailty is the same for recurrent events and the terminal event.

The Multivariate Joint Frailty Model

Multitype recurrent event data arise when one subject experiences two or more different types of recurrent events during the course of the study. As an extension to the joint frailty model for one type of recurrent event and a terminal event (Rondeau et al. 2007), Mazroui and colleagues (2013) introduce joint modeling of two types of recurrent events and a survival outcome. As an example, they estimated the dependency between two different types of recurrent events (locoregional relapses, metastatic relapses) and terminal event (death) after breast cancer. They concluded that the risk of locoregional recurrences was associated with the risk of metastatic recurrences and the risk of metastatic recurrences was also associated with death. As different types of relapses and death were nonindependent events, the hazard of the recurrent events and the hazard of the terminal event (death) needed to be jointly modeled. At a practical level, when predicting the risk of death, accounting for relapses led to better prediction performance. Joint modeling appeared to be suitable for such prediction as well as directly quantifying the correlation between relapses and death.

Following Mazroui et al. (2013), we present the setup and the estimation technique as follows. Let us consider two types of recurrent event times $X_{i j}^{(l)}$ , $j = 1, ..., n_{i}^{(l)}$ for subject $i = 1, ..., N$ ; $l \in 1, 2$ indicates the two types of recurrent events. Let C_i and D_i be censoring and death times. The number $n_{i}^{(l)}$ of observations for recurrent events of type l is a random variable. We denote each individual’s terminal event time as $T_{i}^{*}$ and $T_{i}^{*} = min (C_{i}, D_{i})$ which could be a noninformative censoring C_i or death D_i . $T_{i j}^{(l)} = min (X_{i j}^{(l)}, C_{i}, D_{i})$ corresponds to each follow-up time, $j = 1, ..., n_{i}^{(l)}$ . We consider event times as starting from age $T_{i 0}$ , which is assumed to be 0. $δ_{i j}^{(l)} = I (T_{i j}^{(l)} = X_{i j}^{(l)})$ is a binary indicator for recurrent events, which is 0 if the observation is censored or if the subject had terminal event, and 1 if $X_{i j}^{(l)}$ is observed. Likewise, the death indicator is represented as $δ_{i}^{*} = I (T_{i}^{*} = D_{i})$ . We observe $T_{i j}^{(l)}, T_{i}^{*}, δ_{i j}^{(l)}, δ_{i}^{*}$ .

$N_{i}^{R {(l)}^{*}} (t)$ counts the number of recurrent events of type l for individual i over the interval $(0, t], i = 1, ..., N$ . Due to censoring, we cannot observe the true number of recurrent events experienced by the individual i. Instead, we observe the process $N_{i}^{R (l)} (t) = N_{i}^{R {(l)}^{*}} (min (T_{i}^{*}, t))$ which counts the observed numbers of recurrent events of type l. Similarly, we define the actual and the observed death indicators by time t as $N_{i}^{D^{*}} (t) = I (D_{i} \leq t)$ and $N_{i}^{D} (t) = I (T_{i}^{*} \leq t)$ . Furthermore, let $Y_{i} (t) = I (t \leq T_{i}^{*})$ denote whether or not the individual i is at risk of an event at time t. Over the small interval $[t, t + d t)$ , the number of recurrent events that occur for subject i is $d N_{i}^{R} {(l)}^{*} (t) = N_{i}^{R {(l)}^{*}} ((t + d t)^{-}) - N_{i}^{R {(l)}^{*}} ((t)^{-})$ and the number of observed recurrent events is $d N_{i}^{R (l)} (t) = Y_{i} (t) d N_{i}^{R {(l)}^{*}}$ . Note that $n_{i}^{(l)} = N_{i}^{R (l)} (T_{i}^{*})$ .

The $i^{*} th$ process up to time t is denoted by $H_{i t} = σ {Y_{i} (h), N_{i}^{R (l)} (h), l \in {1, 2}, N_{i}^{D} (h), Z_{i} (h),0 \leq h \leq t}$ , where $Z_{i} (h)$ is a vector of covariates. The intensity processes are jointly dependent through two correlated random effects u_i and v_i , which account for the unobserved heterogeneity, the interrecurrence dependencies, and the dependency between different event types. We assume recurrent events and the terminal event cannot happen at the same time. In addition, we assume that the death precludes the observation of new recurrent events.

The recurrent event intensity processes at time t are, for $l \in (1, 2)$ : $Y_{i} (t) r_{i}^{(l)} (t) d t = P (d N_{i}^{R (l)} (t) = 1 | H_{i t^{-}})$ , where $r_{i}^{(l)} (t) d t = P (d N_{i}^{R {(l)}^{*}} (t) = 1 | Z_{i} (t), u_{i}, v_{i}, D_{i} > t^{-})$ . The death intensity process at time t is $Y_{i} (t) λ_{i} (t) d t = P (d N_{i}^{D} (t) = 1 | H_{i t^{-}})$ , where $λ_{i} (t) d t = P (d N_{i}^{D *} = 1 | Z_{i} (t), u_{i}, v_{i}, D_{i} > t^{-})$ . Finally, we model the intensity functions of counting processes for the two types of recurrent events and the terminal event processes. The multivariate frailty model for two types of recurrent events with a terminal event is given as follows:

{\begin{array}{l} r_{i}^{(1)} (t | u_{i}, v_{i}) = r_{0}^{(1)} (t) e x p (β_{1^{'}} Z_{i 1} (t) + u_{i}) & (recurrence of type 1) \\ r_{i}^{(2)} (t | u_{i}, v_{i}) = r_{0}^{(2)} (t) e x p (β_{2^{'}} Z_{i 2} (t) + v_{i}) & (recurrence of type 2), \\ λ (t | u_{i}, v_{i}) = λ_{0} (t) e x p (β_{3^{'}} Z_{i 3} (t) + α_{1} u_{i} + α_{2} v_{i}) & (terminal event) \end{array}

where $r_{0}^{(l)}, l \in 1, 2$ and $λ_{0} (t)$ are the recurrent and terminal event baseline hazard functions, and $β_{1}, β_{2}, β_{3}$ is the regression coefficient vectors associated with three survival processes. Different covariates, whether time varying or time invariant, could be incorporated into the hazards of the two recurrent events and the terminal event.

The multivariate frailty models are linked together by two correlated Gaussian random effects $u_{i}, v_{i}$ :

(\begin{matrix} u_{i} \\ v_{i} \end{matrix}) \sim N [\begin{matrix} (\begin{matrix} 0 \\ 0 \end{matrix}), & (\begin{matrix} θ & ρ \sqrt{θ η} \\ ρ \sqrt{θ η} & η \end{matrix}) \end{matrix}] .

The random variable u_i is shared by the hazard of recurrent events of type 1 and the hazard of the terminal event. Therefore, the variance of u_i , $θ$ , specifies the variability of the dependencies between occurrences of the recurrent events of type 1 and also interrecurrence dependency. Similarly, the random variable v_i is shared by the hazard of the recurrent events of type 2 and the hazard of the terminal event. The variance of v_i , $η$ , specifies the variability of the dependencies between occurrences of the recurrent events of type 2 and also interrecurrence dependency.

The sign and strength of the parameters $α_{1}$ and $α_{2}$ assess the relationships between two types of recurrent events and the terminal event, respectively. The significance of this parameter estimate informs whether the two types of recurrent events and the terminal event are dependent. Large magnitude of $α_{1}$ ( $α_{2}$ ) illustrates strong dependency between recurrent events of type 1(2) and the terminal event. A high absolute value of the correlation coefficient $ρ$ shows a strong dependency between the two types of recurrent events.

Estimation

In the frailty model, the parameter of interest (regression coefficients, the variance of random variables, and the baseline hazard function) cannot be estimated by maximizing partial likelihood as in Cox models because random variables are added to hazard functions. Several estimation approaches have been proposed: a frequentist approach using the EM algorithm (Klein 1992), a Bayesian approach using Laplace integration (Ducrocq and Casella 1996), or the Markov Chain Monte Carlo algorithm (Clayton 1991). However, these methods are computationally expensive with a large number of iterations. Rondeau, Commenges, and Joly (2003) introduced a semiparametric approach using penalized likelihood.

Unlike parametric models making strong assumptions about the shape of baseline hazard, the semiparametric penalized likelihood methods are very robust. Rondeau et al. (2003) approximate baseline hazard functions by M-splines and estimate parameters $(β_{1}, β_{2}, β_{3}, θ, η, α_{1}, α_{2})$ and the baseline hazard functions $r_{0}^{(t)}, l \in 1, 2$ for the recurrent events or $λ_{0} (t)$ for the terminal event by maximizing the penalized log likelihood. They use minus the converged Hessian of the penalized log likelihood to estimate the variance of the parameters.

Software

Our analyses are implemented in R, with the freely available package frailtypack (Król et al. 2017). This package can be used to fit various joint models for survival events. In particular, multivPenal, a subroutine in frailtypack fits models for two types of recurrent events and a terminal event. The estimates are obtained using a penalized log-likelihood approach, similar to the one developed by Rondeau et al. (2003). The package is available at the Comprehensive R Archive Network (http://cran.r-project.org/package=frailtypack).

The Influence of Employment and Fertility on Retirement Timing

The aging of the labor force in the Western countries stimulates an ever-growing interest among scholars regarding the labor force behavior of older workers. Long-term trends such as soaring rates of unemployment and increases in women’s labor force participation across the life span have implications for the economic behavior of people across the life courses. This interest has recently been extended to include a closer examination of the labor force behavior of older women. Over the second half of the twentieth century, women became far more attached in the workforce than ever before. In Canada, from 1953 to 1990, the labor force participation rate for women grew steadily, rising from about 24 percent in 1953 to 76 percent in 1990. The number reached 82 percent in 2014 (Statistics Canada 2018).

Work and family experiences across the life course, such as employment and childbearing status, have been identified as important determinants influencing women’s labor force decisions in later life and their subsequent retirement timing (Damman, Henkens, and Kalmijn 2015; Finch 2014; Stafford et al. 2019). The work transitions of women in later life are very different from those of men because women have to adjust their work situations to accommodate family responsibilities (Pienta 1999). There are two competing hypotheses in explaining how work–family pathways provide constraints on and opportunities for late-life labor force participation. Attachment hypothesis asserts that paid work participation in earlier adulthood is positively associated with participation in later life (Pienta, Burr, and Mutchler 1994; Stafford et al. 2019). Women who are more attached throughout the life course are more likely to work than women who experienced family-related spells of work interruptions (Finch 2014; Pienta, Burr, and Mutchler 1994).

Alternatively, from an opportunity cost perspective, women with more interruptions due to family caring roles are more likely to work in old age. More time spent in family activities may represent a cost to women’s ability to retire. Women having fewer family responsibilities work more continuously and retire earlier because they build up pension earlier than women having greater family demands (Fisher, Chaffee, and Sonnega 2016; Pienta, Burr, and Mutchler 1994). For some women, time spent raising children comes at the expense of building a stable career. This may mean they attempt to counterbalance these early choices by working more steadily when family demands have decreased (Finch 2014). Our knowledge about the impact of work–family history on retirement timing is mainly from European countries (Damman et al. 2015; Finch 2014; Hank 2004; Stafford et al. 2019). This application seeks to explore these two hypotheses in the Canadian context.

The interlinking of work and family patterns has long been recognized in life-course research (Elder 1994; Elder, Johnson, and Crosnoe 2003). Previous studies have accounted for the links between paid work and family factors in statistical models which only control for summary indicators of employment and parenthood (Hank 2004; Pienta 1999). Few studies combined work and family histories and therefore failed to capture status and timing of transitions in relation to retirement in later life. Joint modeling of several survival processes proves to be a powerful approach to fully utilize retrospective data on work and family status collected across life spans.

Individuals may experience multiple types of recurrent events in their lifetime. For example, both labor force participation and fertility are dynamic processes unfolding over time. Individuals may enter, exit, and reenter the paid labor force at any point in time. Thus, a work trajectory is composed of a sequence of work transitions which may include employment, unemployment, and nonparticipation over adult life. Moreover, work transitions may be linked to changes in martial or parental roles. For example, both number and timing of fertility intentions and behaviors may influence and be influenced by employment at multiple points (Brewster and Rindfuss 2000). We acknowledge the complexity of the transitions from one state to the next and the pathways women follow through the life course as well as the complicated interplay among different life events building over time.

This application aims to examine how employment and childbearing history affect women’s timing of retirement. We contribute to the life course and retirement literature in four respects. First, beyond snapshot information on single events in the life course, we take into account retrospective information on the entire adult life course. Second, prior research has most often analyzed proximate antecedents of retirement and examine work–family pathways during a narrow period of time, we keep track of women’s work and family transitions over a much longer time span from ages 15 to the timing of retirement (or interview time if not retired). Third, we acknowledge the importance of interlinking family processes and employment histories by conceptualizing adult life courses as endogenous survival processes. Fourth, in contrast to the emergence of research of employment and fertility on retirement in European countries and in the United States, we examine this research question in the Canadian context.

Data

We used data from the 2011 Canadian GSS-25 conducted by Statistics Canada. The GSS program collects individual and household-level information on Canadian adults to monitor social changes and people’s well-being (Statistics Canada 2011). The GSS-25 not only records detailed information on retirement decisions, employment, and childbearing histories but also provides standard demographic and socioeconomic data. Thus, it enables an in-depth exploration of how employment and childbearing histories are associated with retirement timing.

In order to capture women with a great degree of labor force volatility, women aged greater than 55 and less than 65 years at the time of interview were selected. This enables us to evaluate the labor force status of women who were approaching the normative age of retirement. We also excluded individuals who never participated in the labor force or worked for less than six months throughout their lives.

As we jointly model two recurrent event processes and one terminal event process, we had three dependent variables (one for each process). The retirement indicator was defined by a woman’s age of retirement. The event indicator for employment was restricted to at most five work spells. It was defined by a woman’s starting age at each of the spells of employment. The event indicator for childbearing was restricted to at most seven children. It was measured by a woman’s age at birth of child 1 to at most child 7 if she had biological child(ren). Exposure times to the risk of the terminal event (retirement) was measured from age 15 until the age of retirement (up to age less than 65). If a respondent did not retire, time was measured from the age 15 to the interview time. Exposure times to the risk of two recurrent events (employment and childbearing) were measured from age 15 until the time of the survey. All subjects were at risk of new recurrent events as long as they were not censored or experienced the terminal event (Cook and Lawless 2007). Once subjects experienced retirement, they were no longer at risk Technical details are given in Supplemental Material (which can be found at http://smr.sagepub.com/supplemental/). Our final study sample included 2,177 female respondents, representing 5,371 employment recurrent events and 3,800 childbearing recurrent events, as well as 430 women undergoing retirement.

As we used retrospective data, we generated life histories beginning at age 15, building in the sequence of employment and childbearing events and followed them until retirement or censoring at the interview time. Time-dependent covariates were subject to change during the life history. For example, if a respondent got married at age 28, then the marital status variable was updated at the time that the marriage began. If a respondent divorced at age 40, the variable was updated again correspondingly. Such covariates can be considered as time dependent and they require special attention to set them up in a counting process style.

The majority of the respondents did not complete their education by age 15. Education is a time-dependent variable which was updated to the highest educational attainment (1 = bachelor’s degree or higher, 2 = some postsecondary education, and 3 = a high school diploma or less) at the date of the respondents’ reported age of education completion.

The GSS-25 collected information on the chronology of marriages. For up to four episodes of marriage, respondents were asked “Age of respondent at start of marriage” as well as “Age of respondent at time of divorce from marriage.” Based on the entry and exit age of marriage, we constructed respondents’ time-dependent marital status variable. Table 1 presents the coding and wording of the dependent and independent variables.

Table 1.

Coding of Variables and Wording of Survey Questions.

Variables	Coding	Wording
Terminal event
Retirement indicator	A woman’s age of retirement	“Age of respondent when she retired.”
Retirement duration	From the age 15 to retirement age (up to age less than 65); from the age 15 to the interview time (if censored)
Recurrent events
Employment indicator	A woman’s starting age at most five spells of employment	“Age of respondent at the beginning of first work period.”“Age of respondent at the beginning of second work period.”“Age of respondent at the beginning of third work period.”“Age of respondent at the beginning of fourth work period.”“Age of respondent at the beginning of fifth work period.”
Employment exposure time	From the age 15 to the interview time
Childbearing indicator	A woman’s age at birth of child 1 to at most child 7 if she had biological child	“Age of respondent at birth of Child_1, if biological child.”“Age of respondent at birth of Child_2, if biological child.”“Age of respondent at birth of Child_3, if biological child.”“Age of respondent at birth of Child_4, if biological child.”“Age of respondent at birth of Child_5, if biological child.”“Age of respondent at birth of Child_6, if biological child.”“Age of respondent at birth of Child_7, if biological child.”
Childbearing exposure time	From the age 15 to the interview time
Time-dependent covariates
Education	Highest level of education:1 = Bachelor’s degree or higher;2 = Some postsecondary education;3 = A high school diploma or less	“Age of respondent at completion of studies.”“Highest level of education obtained by the respondent.”(1 = Doctorate/masters/bachelor’s degree; 2 = Diploma/certificate from community college or trade/technical; 3 = Some university/community college; 4 = High school diploma; 5 = Some secondary/elementary/no schooling)”
Marital status	Married1 = Yes0 = No (divorced or single)	“Age of respondent at start of first marriage.”“Age of respondent at time of divorce from first marriage.”“Age of respondent at start of second marriage.”“Age of respondent at time of divorce from second marriage.”“Age of respondent at start of third marriage.”“Age of respondent at time of divorce from third marriage.”“Age of respondent at start of current marriage.”“Age of respondent at time of divorce from current marriage.”

Table 2 cross-tabulates the number of employment recurrences ranging from 1 to 5 and childbearing recurrences ranging from 0 to 7 for women who retired. For retired women with two children, the number of one, two, and three employment recurrences were 36, 31, and 54, respectively. Among retired women with one employment spell, 58 were childless and 36 had two children.

Table 2.

Cross-tabulation of Number of Children Versus Number of Employment Spells for Retired Women Aged Between 55 and 65.

Number of Children	Number of Employment Spells
Number of Children	1	2	3	4	5	Total
0	58	24	15	7	0	104
1	21	18	15	11	2	67
2	36	31	54	40	15	176
3	8	23	17	11	7	66
4	5	3	1	2	1	12
5	2	1	1	0	0	4
6	1	0	0	0	0	1
7	0	0	0	0	0	0
Total	131	100	103	71	25	430

Note: n = 430.

Results

Survival curve of retirement

As an initial overview of the retirement timing tendencies, we present the survival function of retirement, that is, the probability of retiring at age greater than or equal to $t + 15$ years for t between 0 and 50 years. Of the 2,177 women, 430 retired during the study period and 1,747 were censored. The survival function in Figure 1 decreased rapidly after age 50, corresponding to $t = 35$ . As the respondents in our sample were between 55 and 65 years old, and the oldest retired member observed that the survival curve decreased to 0 just before age 65, corresponding to $t = 50$ .

Figure 1.

Survival curve of retirement.

Recurrent event plot

As an illustration of the recurrent event, we subset four subjects in our sample and displayed their life trajectories (see Figure 2). The figure shows employment and childbearing recurrences before transition to retirement (or censored) for these subjects. Two subjects (numbers 4757, 17648) experienced retirement (represented by a red triangle at the end of the individual’s line segment) and two were censored (numbers 7374, 21594) at the end of the study. Subject number 7374 had the largest number of events. The subject first experienced four employment events and then had two childbearing recurrences and was censored for retirement at the interview time. Subject number 4757 had three employment recurrences at age 19, 32, and 51 and retired at age 54. These examples vividly reflect the complex nature of the multitype recurrent events in that different types of events occur in different orders, and the time between events also differs across subjects.

Figure 2.

Recurrent event plot.

Results with separate models and a joint model

First, we fitted two counting process model of Andersen-Gill, which generalizes the Cox model, for two types of recurrent events (employment and fertility) and fitted a Cox model for terminal event (retirement). Second, we ran a multivariate frailty model with baseline hazard functions approximated by M-splines to model jointly employment recurrences, fertility recurrences, and retirement for women. R codes can be found in Supplemental Material (which can be found at http://smr.sagepub.com/supplemental/). We compared the results from the multivariate frailty model with the results obtained using three separate models. Besides, we were interested in determining the dependencies between two types of recurrent events and the terminal event.

The results are presented in Table 3. The recurrence rate of employment decreased for subjects with less education; the hazard ratio for “Some postsecondary” and “A high school diploma or less” were both below one meaning that employed women with lower levels of education were less likely to be reemployed than respondents obtaining “Bachelor’s degree or higher.” This was consistent with previous research showing the positive impact of education on employment (England, Garcia-Beaulieu, and Ross 2004). There were significant educational differences in the recurrence rate of childbearing. Women with high school diploma or less had higher risk of childbearing than their more educated counterparts. Education was significantly associated with retirement. Women with bachelor’s degree or more were expected to retire earlier than those with lower levels of education.

Table 3.

Separate Models of the Employment and Childbearing Recurrences and Retirement.

Variable	Modalities	Dependent Variable
		Employment Recurrences	Childbearing Recurrences	Retirement
		Coefficient, HR (SD)	Coefficient, HR (SD)	Coefficient, HR (SD)
Education (reference: bachelor’s degree or more)	Some postsecondary	$- 0.196$ , ${0.822}^{* * *}$ (0.044)	$0.521$ , ${1.683}^{* * *}$ (0.060)	$- 0.648$ , ${0.523}^{* * *}$ (0.031)
Education (reference: bachelor’s degree or more)	High school or less	$- 0.454$ , ${0.635}^{* * *}$ (0.043)	$0.556$ , ${1.743}^{* *}$ (0.055)	$- 0.562$ , ${0.570}^{* *}$ (0.085)
Marriage (reference: No)	Yes	$- 0.027$ , $0.973$ (0.031)	$1.892$ , ${6.633}^{* * *}$ (0.052)	$0.353$ , $1.423$ (0.071)
Log likelihood		$- 40, 960.55$	$- 28, 190.79$	$- 3, 076.273$

^a Data are based on the 2011 Canadian General Social Survey.

$* p < 0.05. * * p < .01. * * * p < .001$ (two-tailed test).

Marital status had no significant effect on the recurrence rate of employment. In addition, married women were more likely to have recurrent episodes of childbearing than divorced or single women (hazard ratio [HR] = 6.633). Surprisingly, marital status did not bear a significant association with the hazard of retirement.

The separate models detect only a significant effect of marital status for the recurrence rate of childbearing; it had no significant effect on the recurrence rate of employment or on the hazard of retirement. However, marital status was significantly associated with the recurrence rate of employment, the recurrence rate of childbearing, and the hazard of retirement in the joint model (see Table 4). This illustrated that if the dependencies between these three events of interest were not taken into account, we might omit significant associations. We further observed underestimation of the regression coefficients in the separate models in all three survival processes. For instance, the hazard ratio of marital status was greater in the childbearing process using the joint model compared to the separate model (6.633 vs. 7.056), which did not take into account the terminal event retirement.

Table 4.

Joint Models of the Employment and Childbearing Recurrences and Retirement.

Variable	Modalities	Employment Recurrences	Childbearing Recurrences	Retirement
Variable	Modalities	Coefficient, HR (SD)	Coefficient, HR (SD)	Coefficient, HR (SD)
Education (reference: bachelor’s degree or more)	Some postsecondary	$- 0.170$ , ${0.844}^{* * *}$ (0.044)	$0.524$ , ${1.688}^{* * *}$ (0.060)	$- 0.377$ , ${0.686}^{* * *}$ (0.120)
Education (reference: bachelor’s degree or more)	High school or less	$- 0.357$ , ${0.700}^{* * *}$ (0.042)	$0.556$ , ${1.744}^{* * *}$ (0.054)	$- 0.440$ , ${0.644}^{* * *}$ (0.121)
Marriage (reference: No)	Yes	$0.103$ , ${1.108}^{* * *}$ (0.031)	$1.954$ , ${7.056}^{* * *}$ (0.047)	$0.281$ , ${1.324}^{* *}$ (0.100)
Parameters associated with frailties
		$θ = var (u_{i}) = 0.457$	$η = var (v_{i}) = 0.164$	$α_{1} {= .309}^{* * *}$
		Correlation $ρ = .152$		$α_{2} {= .220}^{* * *}$

^a Data are based on the 2011 Canadian General Social Survey.

$* p < .05. * * p < .01. * * * p < .001$ (two-tailed test).

The added value of the joint model was the frailty parameter estimates (see the lower panel of Table 4). The variance of $u_{i} (θ)$ specified the dependency between occurrences of the recurrent events of type 1 and the terminal event and also the interrecurrence dependency. Similarly, the variance of $v_{i} (η)$ specified the dependency between occurrences of the recurrent events of type 2 and the terminal event and also the interrecurrence dependency and the risk of terminal event. The sign and strength of the dependency between the recurrent event type 1(2) and the terminal event was represented by $α_{1}$ and $α_{2}$ . They were both significantly different from 0, meaning that there were strong dependencies between the risk of employment (fertility) recurrences and the risk of retirement. The parameter estimates $θ$ and $η$ were both not significantly from 0, suggesting that there were no dependencies between the processes explained by the nonobserved factors. The results supported opportunity cost hypothesis suggesting women with more childbearing recurrences and employment recurrences tended to delay retirement because of accumulated savings and benefits (Pienta 1999). Employment histories and childbearing activities over the life course appeared to jointly determine the labor force participation in later life.

The association between the recurrence rate of employment and the recurrence rate of fertility, captured by the correlation coefficient ( $ρ$ ), reflected the interdependence of women’s employment and reproductive histories. By applying joint modeling framework, we demonstrated a reciprocal relationship between the two processes. These frailties parameter estimates highlighted the merits of our proposed model and enable us understand how the intersection of work and family roles impacted the labor force decision of older women.

Conclusions

In this study, we introduced a multivariate frailty model with two correlated random effects to simultaneously model two types of recurrent events with a dependent terminal event. It has focused on the extent to which these processes are interrelated, given that they are subject to joint decision making. Simultaneous modeling of multiple survival processes in joint models offered a number of advantages over separate modeling of each outcome. The model not only explicitly assessed possible dependencies among three survival processes but also corrected for biases introduced into regression estimates due to sources of endogeneity.

We applied this method to examine how the timing of retirement is associated with the risk of employment and risk of childbearing. Rather than focusing on the outcome of static summary of employment years or total number of children ever born, our model considered all episodes of employment or fertility experienced by women from age 15 until the time of retirement (the terminal event) or the time of the survey if retirement had not occurred by the time of the survey. Second, we modeled transitions to employment and fertility and retirement jointly, thus allowing for detecting the endogenous survival processes.

We believe this method will be particularly useful in two research scenarios. The first one is when the focus is on the survival outcome and we wish to account for the effects of recurrent events. For example, health benefits of marriage have long been recognized and extensively studied. However, previous research has yielded inconsistent results for older people (Brown and Wright 2017; H. Liu and Waite 2014). To our knowledge, very little research has considered effects of dynamic and complex marital history with possibilities of several entries to and exits from cohabitations and marriages. Other than counting the number of union transitions, our method can accurately record the timing and sequence of multiple cohabitations and marriages. It will be a powerful tool to answer research questions about how repeatable union transition affects health and longevity in the long run.

The second one is when multiple types of survival data are available. The joint model framework appropriately accounts for potential correlations among them. The method demonstrated in this article provides a more nuanced understanding of correlated event processes, especially the issue of endogeneity in survival analysis. For example, Steele (2011) examined the relationship between employment transitions and births. She found that the number and age of children are associated with the timing of a nonemployed woman’s return to work. Along the same vein, promising research avenues would be exploring relationships between employment transitions and union formation or dissolution.

Future endeavors will be devoted to incorporating techniques for examining other frailty distributions. We have chosen a multivariate normal distribution for the random effects because it is flexible in modeling the covariance structure among different types of event. However, the impact of random effects misspecification in joint model framework warrants further investigation.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Grace Li

Supplemental Material

Supplemental material for this article is available online.

References

Amorim

Leila D. A. F.

Cai

Jianwen

. 2015. “Modeling Recurrent Events: A Tutorial for Analysis in Epidemiology.” International Journal of Epidemiology 44(1):324–33.

Axinn

William G.

Ghimire

Dirgha J.

Smith-Greenaway

Emily

. 2017. “Emotional Variation and Fertility Behavior.” Demography 54(2):437–58.

Balbo

Nicoletta

Barban

Nicola

. 2014. “Does Fertility Behavior Spread among Friends?” American Sociological Review 79(3):412–31.

Blossfeld

Hans-Peter

Hamerle

Alfred

Mayer

Karl Ulrich

. 2014. Event History Analysis: Statistical Theory and Application in the Social Sciences. Hove, England: Psychology Press.

Blossfeld

Hans-Peter

Drobnic

Sonja

, eds. 2001. Careers of Couples in Contemporary Society: From Male Breadwinner to Dual-earner Families. Oxford, England: OUP Oxford.

Box-Steffensmeier

Janet M.

Jones

Bradford S.

. 2004. Event History Modeling: A Guide for Social Scientists. Cambridge, England: Cambridge University Press.

Brewster

K. L.

Rindfuss

R. R.

. 2000. “Fertility and Women’s Employment in Industrialized Nations.” Annual Review of Sociology 26(1):271–96.

Brown

S. L.

Wright

M. R.

. (2017). “Marriage, Cohabitation, and Divorce in Later Life.” Innovation in Aging 1(2):igx015.

Clayton

David G.

1991. A Monte Carlo Method for Bayesian Inference in Frailty Models. Biometrics 47(2):467–485.

10.

Commenges

Daniel

Jacqmin-Gadda

Helene

. 2015. Dynamical Biostatistical Models. Vol. 86. Boca Raton, FL: CRC Press.

11.

Cook

R. J.

Lawless

. 2007. The Statistical Analysis of Recurrent Events. Berlin, Germany: Springer Science & Business Media.

12.

Cox

David R.

1972. “Regression Models and Life-tables (with Discussion).” Journal of Royal Statistical Society 34:187–220.

13.

Damman

Henkens

Kalmijn

. 2015. “Women’s Retirement Intentions and Behavior: The Role of Childbearing and Marital Histories.” European Journal of Population 31(4):339–63.

14.

Dingemans

Mohring

. 2019. “A Life Course Perspective on Working after Retirement: What Role Does the Work History Play?” Advances in Life Course Research 39:23–33.

15.

Ducrocq

Vincent

Casella

. 1996. “A Bayesian Analysis of Mixed Survival Models.” Genetics Selection Evolution 28(6):505.

16.

Elder

G. H.

1994. “Time, Human Agency, and Social Change: Perspectives on the Life Course.” Social Psychology Quarterly 57:4–15.

17.

Elder

G. H.

Johnson

M. K.

Crosnoe

. 2003. “The Emergence and Development of Life Course Theory.” Pp. 3–19 in Handbook of the Life Course, edited by Mortimer

Jeylan T.

Shanahan

Michael J.

. . Boston, MA: Springer.

18.

England

Garcia-Beaulieu

Ross

. 2004. “Women’s Employment among Blacks, Whites, and Three Groups of Latinas: Do More Privileged Women Have Higher Employment?” Gender & Society 18(4):494–509.

19.

Ermisch

John

Steele

Fiona

. 2016. “Fertility Expectations and Residential Mobility in Britain.” Demographic Research 35:1561–84.

20.

Ezell

Michael E.

Land

Kenneth C.

Cohen

Lawrence E.

. 2003. “Modeling Multiple Failure Time Data: A Survey of Variance-corrected Proportional Hazards Models with Empirical Applications to Arrest Data.” Sociological Methodology 33(1):111–67.

21.

Finch

2014. “Why Are Women More Likely Than Men to Extend Paid Work? The Impact of Work-family Life History.” European Journal of Ageing 11(1):31–39.

22.

Fisher

G. G.

Chaffee

D. S.

Sonnega

. 2016. “Retirement Timing: A Review and Recommendations for Future Research.” Work, Aging and Retirement 2(2):230–61.

23.

Hank

2004. “Effects of Early Life Family Events on Women’s Late Life Labour Market Behaviour: An Analysis of the Relationship between Childbearing and Retirement in Western Germany.” European Sociological Review 20(3):189–98.

24.

Hougaard

1995. “Frailty Models for Survival Data.” Lifetime Data Analysis 1(3):255–73.

25.

Klein

John P.

1992. “Semiparametric Estimation of Random Effects Using the Cox Model Based on the EM Algorithm.” Biometrics 48(3):795–806.

26.

Klein

John P.

Moeschberger

Melvin L.

. 2005. Survival Analysis: Techniques for Censored and Truncated Data. Berlin, Germany: Springer Science & Business Media.

27.

Król

Mauguen

Mazroui

Laurent

Michiels

Rondeau

. 2017. “Tutorial in Joint Modeling and Prediction: A Statistical Software for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event.” Journal of Statistical Software 81(1):1–52.

28.

Lillard

Lee A.

1993. “Simultaneous Equations for Hazards.” Journal of Econometrics 56:189–217.

29.

Lillard

Lee A.

Waite

Linda J.

. 1993. “A Joint Model of Marital Childbearing and Marital Disruption.” Demography 30(4):653–81.

30.

Lin

L. A.

Luo

Chen

B. E.

Davis

B. R.

. 2017. “Bayesian Analysis of Multi-type Recurrent Events and Dependent Termination with Non-parametric Covariate Functions.” Statistical Methods in Medical Research 26(6):2869–84.

31.

Liu

Waite

. 2014. “Bad Marriage, Broken Heart? Age and Gender Differences in the Link between Marital Quality and Cardiovascular Risks among Older Adults.” Journal of Health and Social Behavior 55(4):403–23.

32.

Liu

Lei

Wolfe

Robert A.

Huang

Xuelin

. 2004. “Shared Frailty Models for Recurrent Events and a Terminal Event.” Biometrics 60(3):747–56.

33.

Mazroui

Yassin

Mathoulin-Pélissier

Simone

MacGrogan

Gaetan

Brouste

Véronique

Rondeau

Virginie

. 2013. “Multivariate Frailty Models for Two Types of Recurrent Events with a Dependent Terminal Event: Application to Breast Cancer Data.” Biometrical Journal 55(6):866–84.

34.

Myrskylä

Mikko

Fenelon

Andrew

. 2012. “Maternal Age and Offspring Adult Health: Evidence from the Health and Retirement Study.” Demography 49(4):1231–57.

35.

Palloni

Alberto

Beltrán-Sánchez

Hiram

. 2017. “Discrete Barker Frailty and Warped Mortality Dynamics at Older Ages.” Demography 54(2):655–71.

36.

Pessin

Lea

. 2018. “Changing Gender Norms and Marriage Dynamics in the United States.” Journal of Marriage and Family 80(1):25–41.

37.

Petersen

1995. “Models for Interdependent Event-history Data: Specification and Estimation.” Sociological Methodology 25:317–75.

38.

Pienta

A. M.

Burr

J. A.

Mutchler

J. E.

. 1994. “Women’s Labor Force Participation in Later Life: The Effects of Early Work and Family Experiences.” Journal of Gerontology 49(5):S231–39.

39.

Pienta

A. W.

1999. “Early Childbearing Patterns and Women’s Labor Force Behavior in Later Life.” Journal of Women & Aging 11(1):69–84.

40.

Rondeau

Virginie

Commenges

Daniel

Joly

Pierre

. 2003. “Maximum Penalized Likelihood Estimation in a Gamma-frailty Model.” Lifetime Data Analysis 9(2): 139–53.

41.

Rondeau

Virginie

Mathoulin-Pelissier

Simone

Jacqmin-Gadda

Hélène

Brouste

Véronique

Soubeyran

Pierre

. 2007. “Joint Frailty Models for Recurring Events and Death Using Maximum Penalized Likelihood Estimation: Application on Cancer Events.” Biostatistics 8(4):708–21.

42.

Schimmele

Christoph M.

Zheng

. 2016. “Repartnering after Union Dissolution in Later Life.” Journal of Marriage and Family 78(4):1013–31.

43.

Schneider

Daniel

. 2011. “Wealth and the Marital Divide.” American Journal of Sociology 117(2):627–67.

44.

Stafford

Lacey

Murray

Carr

Fleischmann

Stansfeld

McMunn

. 2019. “Work-family Life Course Patterns and Work Participation in Later Life.” European Journal of Ageing 16(1):83–94.

45.

Statistics Canada. 2011. General Social Survey—Family (GSS). Retrieved August 23, 2019 (http://www23.statcan.gc.ca/imdb/p2SV.pl?Function=getSurvey&Id=68196).

46.

Statistics Canada. 2018. The Surge of Women in the Workforce. Retrieved August 23, 2019, (https://www150.statcan.gc.ca/n1/pub/11-630-x/11-630-x2015009-eng.htm).

47.

Steele

Fiona

. 2011. “Multilevel Discrete-time Event History Models with Applications to the Analysis of Recurrent Employment Transitions.” Australian & New Zealand Journal of Statistics 53(1):1–20.

48.

Therneau

T. M.

Grambsch

P. M.

. 2013. Modeling Survival Data: Extending the Cox Model. Berlin, Germany: Springer Science & Business Media.

49.

Vaupel

James W.

Carey

J. R.

Christensen

Johnson

T. E.

Yashin

A. I.

Holm

N. V.

Iachine

I. A.

Kannisto

Khazaeli

A. A.

Liedo

Longo

V. D.

. 1998. “Biodemographic Trajectories of Longevity.” Science 280(5365):855–60.

50.

Vaupel

James W.

Manton

Kenneth G.

Stallard

Eric

. 1979. “The Impact of Heterogeneity in Individual Frailty on the Dynamics of Mortality.” Demography 16(3):439–54.

51.

Wienke

Andreas

. 2010. Frailty Models in Survival Analysis. Boca Raton, FL: CRC Press.

52.

Wrigley-Field

Elizabeth

. 2014. “Mortality Deceleration and Mortality Selection: Three Unexpected Implications of a Simple Model.” Demography 51(1):51–71.