A Joint Model for Longitudinal and Time-to-event Data in Social and Life Course Research: Employment Status and Time to Retirement

Abstract

Longitudinal studies including a time-to-event outcome in social research often use a form of event history analysis to analyse the influence of time-varying endogenous covariates on the time-to-event outcome. Many standard event history models however assume the covariates of interest to be exogenous and inclusion of an endogenous covariate may lead to bias. Although such bias can be dealt with by using joint models for longitudinal and time-to-event outcomes, these types of models are underused in social research. In order to fill this gap in the social science modelling toolkit, we introduce a novel Bayesian joint model in which a multinomial longitudinal outcome is modelled simultaneously with a time-to-event outcome. The methodological novelty of this model is that it concerns a correlated random effects association structure that includes a multinomial longitudinal outcome. We show the use of the joint model on Danish labour market data and compare the joint model to a standard event history model. The joint model has three advantages over a standard survival model. It decreases bias, allows us to explore the relation between exogenous covariates and the longitudinal outcome and can be flexibly extended with multiple time-to-event and longitudinal outcomes.

Introduction

Longitudinal studies with time-to-event outcomes in social research often use a form of event history analysis (EHA) to analyse the influence of (time-varying) covariates on one or several time-to-event outcomes. Applications in social and life course research include a range of event types like divorce (Rosenfeld and Roesler 2019), birth (Baschieri and Hinde 2007), death (Moore and Hayward 1990), job change (Carroll and Mayer 1986; Becker and Blossfeld 2017), interest group survival (Nownes and Lipinski 2005), migrations (Baydar et al. 1990; Henry et al. 2004), and friendship formation and dissolution (Dean et al. 2017).

Both the popular Cox proportional hazards model for the analysis of time-to-event outcomes as well as other types of event history models used to model the effect of time-varying covariates on an event typically assume the relation between the covariate and time-to-event process to be exogenous, meaning that the covariate process does not depend on the survival process (Rizopoulos 2012:44). In other words, the path of an exogenous time-varying covariate up to a specific time is not allowed to be influenced by the occurrence of an event at that time and an anticipated event is thus assumed not to affect the covariate values. Covariates for which this assumption is not met are referred to as endogenous covariates. Failing to appropriately deal with endogenous covariates typically leads to attenuation bias (Ibrahim et al. 2010). Within a larger econometric literature such bias is more commonly referred to as simultaneity bias, that is bias induced by the fact that the time-to-event outcome and time-varying covariate are jointly determined. As an illustration we may consider the effect of anticipated retirement age on late career employment status. Imagine two individuals, A and B, both newly unemployed at age 60. Individual A has paid into an early retirement scheme and knows he/she can retire at age 62 whereas individual B has not paid into such a scheme and can only retire at age 65. One can imagine that their motivation to search for and take up new full time employment differs and therefore the future path of their employment status will most likely differ as well. The likelihood of these individuals taking up new full time employment thus depends on their anticipated retirement age, that is the probability of the event, retirement, occurring at a specific time point. Simply including employment status as a covariate in a Cox model will lead to an underestimation of the effect of employment status on the retirement hazard. A way to deal with covariate endogeneity is to model the time-varying covariate and time-to-event outcome simultaneously as two outcome variables in a joint model for longitudinal and time-to-event data (Lillard and Waite 1993; Faucett and Thomas 1996; Wulfsohn and Tsiatis 1997). Note that in this specific case, with one time-to-event and one categorical longitudinal outcome, multistate models (Cook and Lawless 2018) could also be a viable option for analysis. Apart from covariate endogeneity, joint models also allow for incomplete time-varying covariates as well as measurement error in covariates, two additional problems that cannot be dealt with in the standard Cox model. Joint models have over the past decades received considerable attention in biostatistical research for their ability to deal with measurement error and incomplete longitudinal data as well as their ability to model and quantify endogenous relations between longitudinal covariates and time-to-event outcomes, several longitudinal outcomes and multiple time-to-event outcomes. Since its introduction a multitude of extensions to the standard joint model have been proposed making it a flexible framework for modeling time-to-event and longitudinal variables (see e.g. Hickey et al. (2018) for an extensive overview of available methods for joint models with two or more time-to-event variables and Hickey et al. (2016) for a review of methods for joint models with multiple longitudinal variables).

Whereas in biostatistical and medical research the joint model is well established, in social research joint models have only occasionally been used even though they address an important problem concerning endogenous covariates. In our literature search, we have been able to identify a mere total of six studies in social research that use joint models. In Steele et al. (2013) the relationship between employment transitions and mental health is investigated using a joint model for two longitudinal outcomes. Joint models for a longitudinal and time-to-event outcome are used in Lutz (2014) to explore the relation between occupational activity and the timing of first births and Karimi et al. (2018) for the relation between socio-professional trajectories and cause-specific mortality. Joint models for multiple time-to-event outcomes are used in Lillard and Waite (1993) for the relation between marital conception and marital disruption, Ermisch and Steele (2016) for the relation between fertility expectations and residential mobility and Li et al. (2020) for the relation between employment, fertility and retirement. The majority of social research involving time-to-event outcomes however employs more standard methods that do not allow for covariate endogeneity.

The purpose of the present paper is to extend the social science modeling toolkit beyond the standard event history model to include joint models for longitudinal and time-to-event outcomes. We achieve this goal by first giving a short introduction to joint models for longitudinal and time-to-event outcomes after which we will introduce a novel Bayesian joint model for longitudinal and time-to-event data. This model specifies the relation between the longitudinal and time-to-event outcomes differently from most other (Bayesian) joint models in the literature and can be flexibly extended to include multiple longitudinal and survival outcomes. We will argue that the specific association structure of this joint model may be more suitable in social research, where the longitudinal outcome is of interest in itself and the focus is not solely on the time-to-event outcome. Subsequently, we show the use of the Bayesian joint model, the interpretation of its results and its advantages over a standard event history analysis. We will do so by analysing Danish labour market attachment data in which we jointly model employment and retirement status and comparing the results of this model to those of a standard event history model. We will conclude the paper with a discussion of the advantages and limitations of the joint model as well as some suggestions for further research.

Joint models for Longitudinal and Time-to-Event Outcomes

A joint model for longitudinal and time-to-event outcomes combines at least two submodels, one for the longitudinal and one for the time-to-event outcome. These submodels are referred to as the longitudinal submodel and the survival submodel. The two submodels are linked by means of shared individual-level random effects that can be incorporated in several ways.

The most common approach to joint modeling in the literature is the shared parameter joint model using a current value association structure. In such a model for individuals $i=1,\dots ,N$ and timepoints $t=1,\dots ,n_i$ the longitudinal submodel is specified as a linear mixed-effects model

$$\begin{array}{c}{y}_i(t)={\mathit{x}}_i{(t)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t)}^{\mathrm{\prime }}{\mathit{b}}_i+{ϵ}_i(t)\end{array}$$

(1)

where ${\mathit{x}}_i(t)$ and ${\mathit{z}}_i(t)$ are design matrices $\mathit{\beta }$ is a vector of fixed-effect coefficients, ${\mathit{b}}_i$ contains the individual random effects that are assumed to be normally distributed ${\mathit{b}}_i\sim N(0,{\mathrm{\Sigma }}_{\mathit{b}})$ and ${ϵ}_i(t)$ is a random normal error term with zero mean and variance ${\sigma }^2$ . The survival submodel is specified using a proportional hazards model

$$\begin{array}{c}{\lambda }_i(t)={\lambda }_0(t)\exp ({\mathit{s}}_i{(t)}^{\mathrm{\prime }}\mathit{\gamma }+w_i(t))\end{array}$$

(2)

where ${\mathit{s}}_i(t)$ is a design matrix, $\mathit{\gamma }$ a vector of coefficients, ${\lambda }_0(t)$ is the baseline hazard and $w_i(t)$ is modelled using a current value parameterization $w_i(t)=\alpha ({\mathit{x}}_i{(t)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t)}^{\mathrm{\prime }}{\mathit{b}}_i)$ . The parameter $\alpha$ is the so-called association parameter and quantifies the association between the longitudinal and time-to-event outcome. This association structure is referred to as the current value association structure since it includes an association between the current expected value of the longitudinal outcome and the (log) hazard of the event at time t. This also illustrates one of the differences between the standard Cox model and a joint model. Whereas in a Cox model the observed value of the longitudinal outcome is used as a time-varying predictor, necessitating some form of imputation (most often last observation carried forward), the joint model uses the expected value meaning imputation is not necessary.

Several extensions to this joint model have been proposed. The most common method for incorporating multiple longitudinal outcomes is the specification of a correlated random effects structure (Hickey et al. 2016), other methods have been proposed in Cooper et al. (2007) and Lambert and Vandenhende (2002). In Faucett et al. (1998) and Li et al. (2010) the longitudinal submodel is extended to handle binary or count outcomes. The time-to-event submodel has also been extended in several ways, for example to allow for cause-specific events (Williamson et al. 2008; Huang et al. 2011) or recurring events (Król et al. 2017; Li et al. 2020). The main differences between joint models in the literature however lies in the way they quantify the association between the longitudinal and time-to-event outcome. Apart from the current value association outlined above, lagged value, current and lagged slope, cumulative effects, lagged cumulative effects, weighted cumulative effects, shared random effects, correlated random effects and shared fixed effects associations as well as associations involving interactions with observed data and between longitudinal outcomes are used in the literature (see Table 1 for an overview of some of these association structures). For a more extensive overview of these association structures we refer to Hickey et al. (2016); Hickey et al. (2018) and Brilleman (2018).

Table 1.

Some Association Structures for Joint Models.

Parameterization	Association
Current value	$w_i(t)=\alpha ({\mathit{x}}_i{(t)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t)}^{\mathrm{\prime }}{\mathit{b}}_i)$
Lagged value	$w_i(t)=\alpha ({\mathit{x}}_i{(t-c)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t-c)}^{\mathrm{\prime }}{\mathit{b}}_i)$
Current slope	$\alpha \frac{d^s{\mathit{x}}_i{(t)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t)}^{\mathrm{\prime }}{\mathit{b}}_i}{d{t}^s{\mathit{x}}_i{(t)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t)}^{\mathrm{\prime }}{\mathit{b}}_i}$
Lagged slope	$\alpha \frac{d^s{\mathit{x}}_i{(t-c)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t-c)}^{\mathrm{\prime }}{\mathit{b}}_i}{d{t}^s{\mathit{x}}_i{(t-c)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(t-c)}^{\mathrm{\prime }}{\mathit{b}}_i}$
Cumulative effects	$w_i(t)=\alpha {\int }_0^t{\mathit{x}}_i{(s)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(s)}^{\mathrm{\prime }}{\mathit{b}}_{i}ds$
Weighted cumulative effects	$w_i(t)=\alpha {\int }_0^t\omega (t-s){\mathit{x}}_i{(s)}^{\mathrm{\prime }}\mathit{\beta }+{\mathit{z}}_i{(s)}^{\mathrm{\prime }}{\mathit{b}}_{i}ds$
	where $\omega ()$ is some known weight function
Shared random effects	$w_i(t)=\alpha {\mathit{b}}_i$
Correlated random effects	$w_i(t)=u_i$

The Bayesian model we introduce in the next section makes use of a correlated random effects association structure where $w_i(t)=u_i$ , meaning a frailty parameter $u_i$ is included in the time-to-event submodel. The random effects from the longitudinal submodel and the frailty are then assumed to have a joint multivariate normal distribution:

$$\left[\begin{array}{c}{\mathit{b}}_i\\ u_i\end{array}\right]\sim N\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{\mathrm{\Sigma }}_{\mathit{b}}& {\mathrm{\Sigma }}_{\mathit{b}\mathrm{u}}\\ {{\mathrm{\Sigma }}^{\mathrm{\prime }}}_{\mathit{b}\mathrm{u}}& {\sigma }_u^2\end{array}\right]\right)\cdot$$

The association between the longitudinal and time-to-event submodel is thus captured in the covariances between the frailty and random effects ${\mathrm{\Sigma }}_{\mathrm{bu}}$ .

A Bayesian Correlated Random Effects Joint Model

In this section, we introduce a Bayesian correlated random effects joint model for a categorical longitudinal outcome and one time-to-event outcome. Our choice for a correlated random effects association structure is motivated by its flexibility as well as its better fit to social data types. As exemplified by Elashoff et al. (2008); Huang et al. (2011) and Li et al. (2020) a joint model with a correlated random effects association structure can easily be extended to a competing risk or recurrent events setting. Moreover, by fitting our model using the Stan framework and making code available (see the Estimation section below for more details) researchers wishing to use our model can adapt it to their setting, i.e. including additional outcomes, different prior distributions or models for the longitudinal and time-to-event outcomes, without having to worry about deriving conditional distributions for a Gibbs sampler or having to implement their own Metropolis Hastings algorithms. Finally, the fact that the current value association is the preferred way to implement a joint model in the biostatistical literature, is partly explained by medical research's main interest in the time-to-event outcome, and controlling for a possible influence of exogenous longitudinal variables influencing the time-to-event outcome. In social research however, the longitudinal outcome may be of interest in itself as well. A correlated random effects association structure may in that case be more appropriate, since it does not model a directional effect from the longitudinal outcome onto the time-to-event outcome. In fact, half of the studies in social research employing joint models named in the introduction use a correlated random effects structure.

The methodological novelty of our model is due to the fact that we are to our knowledge the first to use a Bayesian approach to fit a joint model that employs a hierarchical multinomial model for the longitudinal outcome. Such a multinomial structure is useful in social data, where it is often of interest to investigate and compare distinct groups of individuals. This is also exemplified by the six studies employing joint models we named in the introduction. Three of these studies employ a joint model that contains a multinomial longitudinal outcome.

The probability that we observe state $k\in \{1,\dots ,K\}$ in the hierarchical multinomial model for the longitudinal submodel, for individuals $i=1,\dots ,N$ and timepoints $t=1,\dots ,n_i$ equals:

$$\begin{array}{r}{\pi }_{itk}=P(Y_{it}=k)={\displaystyle \frac{1}{1+{\sum }_{h=1}^{K-1}\exp ({\mathit{X}}_{it}\mathit{\beta }+b_{ih})}\mathrm{if}k=K}\\ \begin{array}{c}={\displaystyle \frac{\exp ({\mathit{X}}_{it}\mathit{\beta }+b_{ik})}{1+{\sum }_{h=1}^{K-1}\exp ({\mathit{X}}_{it}\mathit{\beta }+b_{ih})}\mathrm{if}k=1,\dots ,K-1}\end{array}\end{array}$$

(3)

where ${\pi }_{itk}=P(Y_{it}=k)$ is the probability of being in state $k\in \{1,\dots ,K\}$ conditional on the design matrix ${\mathit{X}}_{it}$ , that contains a vector of 1's to be able to estimate a fixed intercept and a set of exogenous covariates s for each individual at each timepoint, and the random intercept $b_{ih}$ for individual $i$ in state $h\in \{1,\dots ,K-1\}$ . The vector $\mathit{\beta }=({\beta }_0,{\beta }_1,\dots ,{\beta }_s)$ is a vector of coefficients for the fixed effects of the exogenous covariates and fixed intercept.

For the time-to-event submodel we use a proportional hazards model with piece-wise constant baseline hazard. In practice, this model can be fit using a hierarchical Poisson regression model for the event indicators $d_{it}$ for individuals $i=1,\dots ,N$ and timepoints $t=1,\dots ,n_i$ (Holford 1980; Laird and Olivier 1981). In case the number of time intervals equals the number of distinct times that events occur, this model is equivalent to the Cox proportional hazard model (Vermunt 1997:96). We fit a hierarchical Poisson log-linear model to model the log pension hazard of the individual i at timepoint $t$ :

$$\begin{array}{c}\log {\mu }_{it}=\log t_{it}+{\mathit{X}}_{it}\mathit{\eta }+\log {\lambda }_t+u_i\end{array}$$

(4)

where ${\mu }_{it}$ is the hazard for individual i at timepoint t, $\log {\lambda }_t$ is the log baseline hazard, $t_{it}$ denotes the time at risk for individual i during the time interval t (the offset), ${\mathit{X}}_{it}$ is a design matrix containing a set of exogenous covariates for each individual at each timepoint, $\mathit{\eta }=({\eta }_1,\dots ,{\eta }_s{)}^{\mathrm{\prime }}$ is a vector of coefficients for the fixed effects of the exogenous covariates and $u_i$ is a frailty parameter (random intercept for each individual). We assume a random walk structure for the baseline hazard as follows: $\log {\lambda }_t={\mu }_{\lambda }+N(\log {\lambda }_{t-1},{\sigma }_{\lambda })$ . Although we could have chosen a more flexible baseline hazard structure, e.g. by specifying an unstructured baseline hazard, the random walk structure adds smoothness to the shape of the baseline hazard. Specifically, the parameter ${\sigma }_{\lambda }$ controls the amount of smoothing with low values representing a high amount of smoothing.

The event history and longitudinal submodels are linked using a joint variance-covariance matrix for the random effects. The random intercepts ${\mathit{b}}_k=b_{1k},\dots ,b_{nk}$ from the longitudinal submodel and the frailties $\mathit{u}=u_1,\dots ,u_n$ from the time-to-event submodel are assumed to follow a multivariate Gaussian distribution $N_{(K-1)+1}(0,\mathrm{\Sigma })$ where

$$\mathrm{\Sigma }=\left[\begin{array}{cc}{\mathrm{\Sigma }}_{\mathit{b}}& {\mathrm{\Sigma }}_{\mathit{b}\mathrm{u}}\\ {{\mathrm{\Sigma }}^{\mathrm{\prime }}}_{\mathit{b}\mathrm{u}}& {\sigma }_u^2\end{array}\right],$$

and where

$${\mathrm{\Sigma }}_{\mathit{b}}=\left[\begin{array}{cc}{\sigma }_{b1}^2& {\rho }_{b1b2}{\sigma }_{b1}{\sigma }_{b2}\\ {\rho }_{b1b2}{\sigma }_{b2}{\sigma }_{b1}& {\sigma }_{b2}^2\end{array}\right],$$

for K = 3 and ${\mathrm{\Sigma }}_{\mathit{b}\mathrm{u}}=({\rho }_{b1u}{\sigma }_{b1}{\sigma }_u,{\rho }_{b2u}{\sigma }_{b2}{\sigma }_u{)}^{\mathrm{\prime }}$ . This means that we can assess the association between the longitudinal and event history process using the covariances in ${\mathrm{\Sigma }}_{\mathit{b}\mathrm{u}}$ . E.g.: if ${\rho }_{b1u}{\sigma }_{b1}{\sigma }_u>0$ this means larger values of $b_{i1}$ tend to coincide with larger values of $u_i$ and if ${\rho }_{b1u}{\sigma }_{b1}{\sigma }_u<0$ this means larger values of $b_{i1}$ tend to coincide with smaller values of $u_i$ and vice versa. In the context of the joint model a positive covariance should be interpreted as: a larger individual probability of being in a certain state k coincides with a larger individual hazard.

Estimation

The Bayesian joint model is estimated using a Hamiltonian Monte Carlo (HMC) sampler (Neal 2011; Betancourt and Girolami 2015) specified in the Stan language (Carpenter et al. 2017) in R (Stan Development Team 2019) (code is included in the Appendix). We specify independent uninformative $N(0,{\sigma }_0=10000)$ priors for the fixed effects, ${\mathit{\beta }}_k$ and $\mathit{\eta }$ . The prior for $\log {\lambda }_{t=1}$ and hyperpriors for ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }$ are specified as $N(0,1)$ priors and the prior for $\log {\lambda }_{t>1}$ equals $N(\log {\lambda }_{t-1},{\sigma }_{\lambda })$ . We assume $N(0,1)$ hyperpriors for the means of the random effects and LKJ (Lewandowski Kurowicka Joe) hyperpriors (Barnard et al. 2000) for their correlation matrix. The code to implement this model is included in the Appendix.

Analysis of the Labour Market Attachment Data

In this section we introduce the labour market attachment data and analyse it using our Bayesian joint model. The research question of interest is how a reform in early retirement pension (ERP) affects occupational life-trajectories, especially regarding the relation between employment status and retirement timing.

We will interpret the results from the joint model model and compare these to results from a standard event history model. The standard event history model is equivalent to the time-to-event submodel of the joint model specified in (4). The fit of the standard event history model and the joint model will be evaluated on the complete labour market attachment data by comparing the predicted survival curves of both models to the proportions of retired individuals in the complete data.

Labour market Attachment Data

In Denmark there are three main state-funded retirement schemes that pay monthly benefits to those eligible: old age pension (OAP), early retirement pension and disability pension (DP). Whereas old age pension is available to those over 65, employees that meet certain conditions ¹ can opt to make use of the early retirement pension. In 2006 the Danish government announced to gradually alter the age at which the early retirement pension would be available. The reform means that for individuals born in 1953 or earlier, the early retirement age lies at 60 whereas for individuals born in 1954 or later this age is gradually increased from 60.5 for individuals born between 01-01-1954 and 01-07-1954 to computed in relation to life expectancy for individuals born in 1963 and later.

The labour market attachment data contains individual data on monthly labour market status as well as sex and education level for three cohorts, two unaffected (born in 1953 and 1950) and one affected (born in 1954) by the ERP reform, obtained from the Danish population registers. We follow the individuals from these cohorts from the month they turn 58 (in 2008, 2011 and 2012 respectively) up to the month they turn 63 (in 2013, 2016 and 2017 respectively) to show the effect of the reform on the occupational life-trajectories over time. This leads to a dataset with 212,655 individuals (71,459, 71,245 and 69,951 in the 1950, 1953 and 1954 cohort respectively) measured over 60 months. We will refer to this data as the labour market attachment data. Labour market status is divided into four categories: ‘full-time employed’, ‘part-time employed’ and ‘outside the labour market’ refer to the employment status whereas ‘receiving (any type of) pension’ represents the retirement status (see the Appendix for further details on how these categories were created). Figure 1 shows the full labour market attachement data (including deceased individuals). This means that we can follow the transition of Danish individuals from being employed to retirement, and construct what we call occupational life-trajectories. To reduce computation time of the HMC sampler, we take a stratified (on sex and education level) sample from the larger labour market attachment data leading to subsamples of 610, 609 and 596 individuals from the 1950, 1953 and 1954 cohort respectively. On a Windows server with an Intel Xeon CPU E5-2660 v3 processor with 768 GB RAM the analyses took approximately 4.5 h to finish.

Figure 1.

Attachment to the labour market over time for three cohorts. The y-axis shows the proportion of the population with a particular attachment (‘fulltime employment’, ‘parttime employment’, ‘outside the labour market’ or ‘receiving pension’). In addition to labour market attachment the proportion of individuals that have died is also shown.

In pension reforms such as the reform in ERP one of the main goals is to increase the employment rate of older individuals in order to alleviate the decline of the working-age population and take into account increasing life expectancy. Increasing the employment rate in this way will have a positive fiscal effect by extending pension contribution periods and reducing the number of individuals receiving pension simultaneously. However, increasing the pension age may have a negative impact if older individuals are not able to work longer and instead end up receiving unemployment benefits. The research question of interest with regard to the labour market attachment data is therefore how the reform in ERP affects occupational life-trajectories, especially regarding the relation between employment status and retirement timing. Additionally, we investigate the effect of education level and gender on the occupational life-trajectories.

Joint Model Results

With regard to the labour market attachment data we are interested in the effect of a policy change in ERP on the occupational life-trajectories. To investigate this effect we split up the attachment in two distinct variables, employment status (longitudinal outcome) and retirement status (time-to-event outcome) that we included as outcome variables in three separate Bayesian joint models, each fit to a subset of one of the three different cohorts of the Danish population. In all models sex (male = 1) and education level (post-secondary or tertiary education = 1) are used as binary covariates. The vectors $\mathit{\beta }$ and $\mathit{\eta }$ from (3) and (4) thus consists of three (two covariates and intercept) and two elements respectively. In each of our Bayesian joint models we set up our Stan sampler to use four chains of 4000 (including 1000 burn-in) iterations each and set a lag of 4 to prevent possible autocorrelation.

The joint model allows us to investigate the effect of the policy change in ERP in several ways. Firstly, we can investigate differences between the three cohorts in the effect of sex and education level on retirement timing. Table 2 shows posterior summaries of the hazard ratios of the exogenous covariates in the time-to-event submodel for the three cohorts. We see that females are about 1.5 times more likely to retire than males in the 1950 and 1953 cohort. In the 1954 cohort however, females are estimated to retire 1.9 times more likely than males. Furthermore, individuals with a post-secondary or tertiary education are 0.64 times as likely to retire as individuals with only primary or secondary education in the 1954 cohort whereas in the other two cohorts the likelihood to retire is not affected by education level.

Table 2.

Posterior Mean and 95% Credible Intervals (CI) of the Hazard Ratios, exp(η), of the Exogenous Covariates (Sex and Education Level) in the Time-to-Event Submodel of the Joint Model for the Three Cohorts. The First Line (Sex) Represents Females vs. Males and the Second Line (Education) Represents Post-Secondary and Tertiary vs. Primary and Secondary Education.

	1950 cohort		1953 cohort		1954 cohort
	Mean	CI	Mean	CI	Mean	CI
Sex	1.47	(1.15, 1.85)	1.52	(1.18, 1.93)	1.91	(1.45, 2.50)
Education	0.78	(0.59, 1.00)	0.78	(0.59, 1.01)	0.64	(0.47, 0.85)

Secondly, we can investigate differences between the three cohorts with regard to the probability of having a certain employment status and the effects of sex and education level on this status. Table 3 shows posterior summaries for the log odds of being employed (fulltime or parttime) versus being outside the labour market. In all three cohorts the odds of being fulltime employed vs. outside the labour market is lower for females than for males. The odds of being parttime employed vs. outside the labour market are higher for females than for males in the 1953 cohort. The odds of being fulltime or parttime employed vs. being outside the labour market are higher for individuals with a high education except for parttime employed individuals in the 1950 cohort. Table 4 shows the probability of having a certain employment state. The probability of being fulltime employed increases over time for all groups except males with a high education level. Whereas the probability of fulltime employment increases from the 1950 to the 1954 cohort, the probability of parttime employment decreases for all groups except males with a higher education (see rows 5 to 8 in Table 4). Note that the increases/decreases in time are not linear for all sex and education subgroups e.g. for females with a high education level the probability of parttime employment increases from the 1950 to 1953 cohort and decreases again from the 1953 to 1954 cohort.

Table 3.

Posterior Mean and 95% Credible Intervals (CI) of the Intercept ${\beta }_0$ and the Coefficients for the Effects of Sex and Education Level ${\beta }_1$ and ${\beta }_2$ in the Joint Model for the Three Cohorts. The Value $\exp (\beta )$ Gives the Corresponding Odds Ratios. The Reference Categories for Sex and Education are Males and Primary and Secondary Education Respectively.

		1950 cohort			1953 cohort			1954 cohort
		Mean			Mean			Mean
		$\beta$	$\exp (\beta )$	CI	$\beta$	$\exp (\beta )$	CI	$\beta$	$\exp (\beta )$	CI
Fulltime employment	intercept	2.65	14.15	(1.94, 3.38)	2.11	8.24	(1.38, 2.95)	3.02	20.41	(2.14, 3.89)
	sex	−2.23	0.11	(−3.27, −1.22)	−1.79	0.17	(−2.80, −0.68)	−2.07	0.13	(−3.20, −0.93)
	education	1.83	6.24	(0.71, 2.97)	1.98	7.24	(0.81, 3.20)	2.29	9.83	(1.03, 3.50)
Parttime employment	intercept	1.51	4.53	(0.99, 2.06)	0.60	1.82	(−0.13, 1.28)	1.20	3.31	(0.48, 1.99)
	sex	0.66	1.93	(−0.07, 1.35)	0.95	2.60	(0.10, 1.85)	0.30	1.35	(−0.65, 1.22)
	education	0.54	1.72	(−0.26, 1.39)	1.52	4.56	(0.50, 2.57)	1.75	5.75	(0.68, 2.77)

Table 4.

Posterior Mean and 95% Credible Intervals (CI) of the Probabilities of Having a Certain Employment status Conditional on Sex (Male/Female) and Education Level (low/high) in the Joint Model for the Three Cohorts.

		1950 cohort		1953 cohort		1954 cohort
		Mean	CI	Mean	CI	Mean	CI
Fulltime	male + low	0.71	(0.59, 0.82)	0.74	(0.62, 0.84)	0.82	(0.71, 0.90)
Employment	male + high	0.90	(0.81, 0.96)	0.86	(0.73, 0.94)	0.90	(0.80, 0.96)
	female + low	0.14	(0.08, 0.23)	0.20	(0.11, 0.32)	0.32	(0.18, 0.48)
	female + high	0.37	(0.20, 0.57)	0.31	(0.15, 0.49)	0.49	(0.29, 0.69)
Parttime	male + low	0.23	(0.14, 0.34)	0.17	(0.10, 0.26)	0.14	(0.08, 0.23)
Employment	male + high	0.09	(0.04, 0.17)	0.13	(0.05, 0.24)	0.09	(0.04, 0.19)
	female + low	0.77	(0.67, 0.85)	0.66	(0.54, 0.78)	0.55	(0.40, 0.69)
	female + high	0.58	(0.40, 0.74)	0.65	(0.47, 0.81)	0.49	(0.30, 0.68)
Outside the	male + low	0.05	(0.03, 0.09)	0.09	(0.07, 0.16)	0.04	(0.02, 0.09)
Labour	male + high	0.01	(0.00, 0.03)	0.02	(0.00, 0.04)	0.01	(0.00, 0.01)
Market	female + low	0.09	(0.05, 0.15)	0.15	(0.07, 0.24)	0.13	(0.06, 0.23)
	female + high	0.04	(0.02, 0.08)	0.03	(0.01, 0.07)	0.02	(0.01, 0.05)

Finally, the joint model allows us to explore the relation between employment status and retirement timing for the three cohorts. The positive correlation ${\rho }_{b1b2}$ in the first row of Table 5 indicates that individuals that are more likely to be fulltime employed are also more likely to be parttime employed and vice versa. The negative correlations ${\rho }_{b1u}$ and ${\rho }_{b2u}$ in the second and third rows of Table 5 indicate that individuals that are more likely to be employed have a lower pension hazard and vice versa. Individuals under employment are thus less likely to transition into retirement. Individuals under parttime employment are however more likely to transition into retirement than individuals under fulltime employment, reflected in the weaker negative correlations in Table 5. The three last rows of Table 5 contain the standard deviations of the random effects and frailty for the three cohorts. The standard deviations indicate that individuals became more alike in their pension hazard over time, from the 1950 to the 1954 cohort ( ${\sigma }_u$ decreases), but less alike in their employment trajectories ( ${\sigma }_{b1}$ and ${\sigma }_{b2}$ increase).

Table 5.

Posterior Mean and 95% Credible Intervals (CI) of the Random Effect Correlations and Standard Deviations in the Joint Model for the Three Cohorts. The Parameter $b_1$ is the Random Intercept for Fulltime Employment and $b_2$ is the Random Intercept for Parttime Employment.

	1950 cohort		1953 cohort		1954 cohort
	Mean	CI	Mean	CI	Mean	CI
${\rho }_{b1b2}$	0.61	(0.53, 0.69)	0.64	(0.56, 0.71)	0.65	(0.56, 0.72)
${\rho }_{b1u}$	−0.79	(−0.95, −0.59)	−0.82	(−0.96, −0.62)	−0.77	(−0.95, −0.52)
${\rho }_{b2u}$	−0.36	(−0.59, −0.11)	−0.53	(−0.74, −0.32)	−0.50	(−0.77, −0.20)
${\sigma }_{b1}$	5.68	(5.17, 6.24)	5.97	(5.43, 6.56)	6.32	(5.69, 6.98)
${\sigma }_{b2}$	3.90	(3.55, 4.29)	4.90	(4.45, 5.42)	5.21	(4.70, 5.75)
${\sigma }_u$	0.50	(0.32, 0.70)	0.57	(0.39, 0.79)	0.44	(0.26, 0.64)

Standard Event History Results

To investigate the effect of sex and education level on retirement timing as well as the association between labour market status on retirement timing in the standard event history model we look at the hazard ratios in Table 6. The hazard ratio for sex only indicates an effect in the 1953 (marginal) and 1954 cohort and the hazard ratio for education level indicates a marginal effect in the 1954 cohort. Furthermore, the hazard ratios for the employment status indicate that the pension hazard is lower for individuals that are fulltime or parttime employed in all three cohorts. We cannot investigate differences between the three cohorts with regard to the probability of having a certain employment status since we use employment status as an exogenous predictor instead of as an endogenous outcome variable in the standard even history model.

Table 6.

Posterior Mean and 95% Credible Intervals (CI) of the Hazard Ratios, exp(η), of the Exogenous Covariates (Sex, Education Level and Employment status (Fulltime and Parttime Employment)) in the ‘Standard’ Event History Model for the Three Cohorts.

	1950 cohort		1953 cohort		1954 cohort
	Mean	CI	Mean	CI	Mean	CI
Sex	1.21	(0.95, 1.54)	1.30	(1.01, 1.65)	1.65	(1.25, 2.15)
Education	1.01	(0.74, 1.31)	0.97	(0.74, 1.25)	0.75	(0.55, 0.99)
Fulltime employment	0.13	(0.08, 0.18)	0.20	(0.14, 0.27)	0.20	(0.14, 0.28)
Parttime employment	0.47	(0.35, 0.61)	0.54	(0.41, 0.71)	0.58	(0.43, 0.77)

Model comparison

The results from the standard event history model are slightly different from those of the joint model. The joint model seems to have slightly more power to pick up the effects of sex and education on retirement timing (see the confidence intervals of the hazard ratios in Tables 2 and 6). Second, over time, from the 1950 to the 1954 cohort, the hazard ratios of fulltime and parttime employment in the standard model increase meaning that employed individuals have a higher pension hazard in 1954 than in 1950. This conclusion is opposite to that of the joint model, in which the negative correlations between frailty and random intercept either increase or stay at the same level, indicating the same or a stronger negative relation between employment and retirement. This suggests that not taking the longitudinal data into account can lead to different conclusions. Overall however, the direction of the hazard ratios (and correlation coefficients) of sex, education and employment status, with females, individuals with a shorter education and unemployed individuals being more likely to retire is the same in both the standard and joint model.

Finally, to evaluate the fit of the models we compute Brier scores (BS) for both the joint model and standard event history model. The Brier score can be used to evaluate the predictive accuracy of a survival function at each individual timepoint. It ranges from zero to one, where a zero value represents perfect fit. In a survival setting the Brier score is computed as:

$$BS(t)={\displaystyle \frac{1}{N}\sum _{i=1}^N(1_{T_i>t}-\hat{S}(t,{\overline{x}}_i){)}^2}$$

where N Is the sample size, $T_i$ is the event time for individual i and $\hat{S}(t,{\overline{x}}_i)$ is the predicted survival function at timepoint $t.$ A version that accounts for right-censoring is given in Graf et al. (1999) and Gerds and Schumacher (2006). Computed Brier scores based on the individual estimated posterior mean survival scores as well as the lower and upper bound of the 95% credible interval of the individual survival scores are shown in Figure 2. Even though in theory estimates from the ‘standard’ model may be biased, we see that in case of the labour market attachment data, Brier scores do not differ a lot.

Figure 2.

Brier scores for the joint model and standard event history model (lower values are better). The solid lines are the Brier scores computed using the estimated individual posterior mean survival scores. The upper and lower bounds of the shaded region are the Brier scores for the joint and standard event history model computed using the upper and lower bounds of the 95% credible interval of the individual survival scores.

In order to compare the standard event history model to the joint model in more detail we conduct a small simulation study in the next section.

Simulation study

In order to show the endogeneity bias when fitting standard event history models as well as the effects of censoring on the model estimates. Table 7 shows an overview of the parameter values and sample sizes for the seven designs that were compared in the simulation. In each design, 100 datasets were simulated using a correlated random effects joint model in which the survival outcome is a hierarchical poisson log-linear model with random walk baseline and a hierarchical logistic model is used for the longitudinal outcome. We use the same binomial covariates with success probabilities of 0.47 and 0.3 respectively for both outcomes. Code for the simulation study can be found in the Online Supplement.

Table 7.

Simulation Designs.

Design	$N$	$n_i$	${\eta }_1$	${\eta }_2$	${\beta }_0$	${\beta }_1$	${\beta }_2$	${\mu }_{\lambda }$	${\sigma }_{\lambda }$	${\rho }_{bu}$	Censoring
1	300	30	0.5	−0.3	1.5	−0.5	2	−1	0.45	−0.8	0%	–
2	600	30	0.5	−0.3	1.5	−0.5	2	−1	0.45	−0.8	0%	–
3	1200	30	0.5	−0.3	1.5	−0.5	2	−1	0.45	−0.8	0%	–
4	600	30	0.5	−0.3	1.5	−0.5	2	−1	0.45	−0.8	5%	independent
5	600	30	0.5	−0.3	1.5	−0.5	2	−1	0.45	−0.8	10%	independent
6	600	30	0.5	−0.3	1.5	−0.5	2	−1	0.45	−0.8	15%	independent
7	600	30	0.5	−0.3	1.5	−0.5	2	−1	0.45	−0.8	5%	dependent

Table 8 shows the bias and standard deviation of the posterior means of the estimates. The bias is computed as the average deviation of the posterior mean estimate from the real value. Positive and negative values indicate over- and underestimation, respectively. Note that for each design we fit three models: a joint model, a standard event history model, and a standard event history model that includes the binary longitudinal outcome as a time-varying covariate. From Table 8 we conclude that there is a clear bias in the etimates of the standard event history models in all settings. The joint model has least bias as expected. Furthermore, the bias decreases with increasing sample size. Censoring increases bias in ${\eta }_1$ and ${\eta }_2$ . Surprisingly, censoring does not seem to negatively affect the estimates of ${\mu }_{\lambda }$ , ${\sigma }_{\lambda }$ and ${\rho }_{bu}$ . Estimates for ${\mu }_{\lambda }$ actually seem to improve in models with censoring.

Table 8.

Bias and Standard Deviation of the Estimated Posterior Means.

Analysis	Design	${\eta }_1$	${\eta }_2$	${\beta }_0$	${\beta }_1$	${\beta }_2$	${\mu }_{\lambda }$	${\sigma }_{\lambda }$	${\rho }_{bu}$
Joint	1	−0.17 (0.61)	−0.23 (0.69)	0.05 (0.13)	0.03 (0.13)	0.04 (0.22)	0.35 (0.75)	−0.23 (0.10)	0.15 (0.14)
	2	−0.04 (0.47)	−0.04 (0.58)	0.05 (0.09)	−0.02 (0.10)	0.03 (0.16)	0.32 (0.90)	−0.12 (0.11)	0.08 (0.09)
	3	0.00 (0.30)	0.01 (0.40)	0.02 (0.06)	−0.01 (0.07)	0.00 (0.10)	0.28 (1.08)	−0.08 (0.10)	0.02 (0.07)
	4	−0.14 (0.47)	−0.18 (0.57)	0.05 (0.09)	0.01 (0.10)	0.03 (0.15)	0.08 (0.93)	−0.15 (0.10)	0.08 (0.10)
	5	−0.14 (0.47)	−0.12 (0.63)	0.05 (0.09)	0.01 (0.09)	0.03 (0.17)	−0.02 (0.89)	−0.15 (0.10)	0.07 (0.08)
	6	−0.16 (0.49)	−0.17 (0.51)	0.07 (0.10)	0.00 (0.11)	0.03 (0.16)	−0.20 (0.85)	−0.15 (0.11)	0.07 (0.09)
	7	−0.08 (0.49)	0.13 (0.55)	0.04 (0.11)	−0.00 (0.09)	0.01 (0.18)	0.11 (0.83)	−0.14 (0.10)	−0.06 (0.08)
Standard	1	0.52 (0.18)	1.32 (0.18)	–	–	–	2.02 (0.19)	0.57 (0.18)	–
	2	0.54 (0.18)	1.34 (0.17)	–	–	–	2.04 (0.18)	0.59 (0.18)	–
	3	0.51 (0.17)	1.31 (0.17)	–	–	–	2.01 (0.17)	0.56 (0.17)	–
	4	0.48 (0.17)	1.28 (0.17)	–	–	–	1.98 (0.17)	0.53 (0.17)	–
	5	0.46 (0.20)	1.26 (0.20)	–	–	–	1.96 (0.20)	0.52 (0.20)	–
	6	0.43 (0.16)	1.23 (0.16)	–	–	–	1.93 (0.16)	0.48 (0.16)	–
	7	0.55 (0.16)	1.35 (0.16)	–	–	–	2.05 (0.16)	0.60 (0.16)	–
Standard time-varying	1	0.59 (0.19)	1.39 (0.19)	–	–	–	2.10 (0.19)	0.64 (0.19)	–
	2	0.60 (0.18)	1.40 (0.18)	–	–	–	2.10 (0.18)	0.65 (0.18)	–
	3	0.58 (0.16)	1.37 (0.16)	–	–	–	2.08 (0.16)	0.63 (0.16)	–
	4	0.54 (0.17)	1.34 (0.17)	–	–	–	2.04 (0.17)	0.59 (0.17)	–
	5	0.54 (0.20)	1.34 (0.20)	–	–	–	2.03 (0.20)	0.59 (0.20)	–
	6	0.52 (0.16)	1.32 (0.16)	–	–	–	2.02 (0.16)	0.57 (0.16)	–
	7	0.61 (0.17)	1.41 (0.17)	–	–	–	2.11 (0.17)	0.66 (0.17)	–

Discussion

In this paper we have expanded the social science modelling toolkit by introducing a Bayesian joint model for longitudinal and time-to-event data and illustrating the advantages of this model over standard event history analysis by analyzing the labour market attachment data.

As expected, the comparison between the joint model and a standard event history analysis shows that the joint model better fits the labour market attachment data. Additionally, there are some slight differences regarding the effect of employment on the pension hazard. Even though the application of the joint model in the labour market data does not seem to lead to very strong differences in conclusions, this may not be the case for all data. As argued in the introduction, the main reason to prefer the joint model over standard methods, is the fact that it leads to less statistical concerns (about bias) in case of endogenous time-varying covariates than standard event history models. A secondary reason to prefer the joint model, concerns the fact that we are able to include predictors of the time-varying covariate into the model. This is especially advantageous in social research, where the interest may often not solely lie with the time-to-event outcome but where there may be an interest in the time-varying covariate in itself as well.

In some recent studies on retirement timing instead of using time-to-event analyses, regression discontinuity models are used, e.g. in Geyer and Welteke (2019) and Nielsen (2019). In such models however the focus lies on variables that measure retirement at/before a specific age, in which case censoring is not problematic in the sense that this is an observed measure for all individuals. Using such a model in our case may cause problems due to censoring, e.g. the estimated pension hazards may be biased. Also, as mentioned in the introduction, multistate models may be a viable alternative to the joint model in case of the labour market data used in this paper. However, this is only true because the labour market data contains a time-to-event outcome in combination with a categorical longitudinal outcome. If we decide to extend the model with an additional continuous longitudinal outcome, say montly income, joint models can incorporate such outcomes without having to categorize them as multistate models would.

Although the model we introduce in this paper has a major advantage over standard EHA when it comes to the inclusion of endogenous covariates it also comes at a price. Firstly, the computation speed of the joint model. Even though the number of individuals in the subsample of the data we analysed is not very large and we only need to estimate two random effects and one frailty, the computation speed is quite slow. This is partly due to the long observation period, the fact that we use an iterative Bayesian method to estimate the model and the fact that estimation for multinomial distributions has not been optimized in Stan. There are, however, several possibilities to increase the computation speed of the model. Firstly, instead of using iterative Bayesian methods such as HMC sampling we may instead resort to approximate Bayesian methods, e.g. we could use the variational Bayesian methods for approximating the posterior instead of sampling from it that are implemented in Stan (Kucukelbir et al. 2015). Secondly, we can resort to methods for parallel Bayesian computation, e.g. using Consensus Monte Carlo (Neiswanger et al. 2014) or Wasserstein Posterior (Srivastava et al. 2018) methods, for realizing a speed up in computation time. Lastly, we may reparameterize the model or use different methods for specifying the submodels that lead to faster computations (see Bradley et al. (2018) for an example).

Another caveat concerns the fact that although the random effect in the longitudinal submodel describes the employment trajectory of each individual, by modeling the probability of being fulltime-employed, parttime employed or outside of the labour market over the entire trajectory, this summary cannot compete with the more detailed typification that could have been obtained through so-called ‘holistic’ approaches, e.g. sequence analysis or (latent) growth models (Billari and Piccarreta 2005; Piccarreta and Studer 2019). Sequence analysis methods would for example be able to include both the duration and pattern of different employment trajectories as factors influencing different types of employment trajectories. The addition of lagged covariates to the joint model may allow us to bring the description of the employment trajectory of an individual to a level more close to that of sequence analysis methods. Note however that in sequence analysis it is much more complicated to assess the effect of additional (exogenous) covariates.

The joint model introduced in this paper only considers one time-to-event outcome and one longitudinal outcome. We have for example not included death as a competing risks outcome. This can, however, be dealt with in a relatively simple way, as the joint model is flexible and it is straightforward to add additional longitudinal and/or time-to-event outcomes due to its implementation in Stan. To include additional longitudinal outcomes we simply add submodels for additional longitudinal outcomes and include the random effects from these submodels to the variance-covariance matrix. The complexity of adding time-to-event outcomes depends on the association structure we assume. If we assume simultaneous hazards this issue can be solved in the same way as for adding additional longitudinal outcomes. If we assume competing risks we can add submodels for the additional time-to-event outcomes as well but in that case we need to take care when interpreting the parameter estimates (Andersen et al. 2012). Note that the joint model presented in this paper allows for censoring, e.g. due to emigration. Treating death as censoring leads to a model in which the event (early retirement) is allowed to take place sometime after death. Despite the fact that such a potential outcome space after death may be of interest in some applications, the results from such models are difficult to interpret as they reflect a situation that is impossible in the real world.

In conclusion, we have both introduced and shown the use of a Bayesian joint model for longitudinal and time-to-event data. The joint model has clear advantages over a standard survival model. It decreases bias in estimation, it allows us to explore the relation between exogenous covariates and the time-to-event as well as the longitudinal outcome and it allows for flexible extensions in terms of including further longitudinal and time-to-event variables in the same model.