Bridging the Gap: Introducing Joint Models for Longitudinal and Time-to-Event Data in the Social Sciences

An Illustrative Example: Time-to-Event Model and Joint Model in Comparison

In the example we are interested in — the risk of marriage dissolution (event) and how covariates influence this risk — we may include covariates that are changing over time such as subjective marital satisfaction (TVC). Research interest lies in the influence of the trajectory of marital satisfaction and the risk of marriage dissolution. Assume marital satisfaction has been captured multiple times over the years and information on the start and possible end of the marriage are available. Figure 1 serves as an illustration for the fictional example using one single individual: the upper trajectory depicts his/her marital satisfaction measure, where the points correspond to the measurements in time. The measurement points are connected via the smooth trajectory function (solid line). The lower part represents the estimated hazard rate for marriage dissolution over time. In this example higher values of satisfaction (upper panel) go hand in hand with smaller estimated hazard rates, i.e., lower risk of ending the relationship (lower panel).

Click to enlarge

Figure 1

Scheme of Two Related Processes and Possible Modelling Strategies

Note. Upper panel: Measurements of the time-varying covariate (longitudinal process). Lower panel: Risk of having an event (time-to-event process). In a joint model the influence of the trajectory in the upper panel on the risk is estimated by an association parameter α and both processes can be modelled as functions of (shared) covariates.

In order to include the time-varying covariate marital satisfaction into a time-to-event model, one could use the classical TVC approach. This approach yields several problems since it assumes that the respective variable, (1) does not change between the observation times and, (2) is exogenous. These two assumptions are particularly questionable for volatile and self-reported variables. The first assumption would result in a step-function (dashed line in upper panel of Figure 1) between the measurements in the upper panel instead of the smooth path. Especially for infrequently measured variables with long periods between two measures this approach may model the trajectory inappropriately.

Additionally, as marital satisfaction does neither evolve from a stochastic process, which is external to the individual under study (ancillary TVC), nor can be calculated as a defined covariate, it can be called an endogenous TVC and therefore violates the second assumption. Figure 2 may be an indicator for state dependence of marital satisfaction showing the smoothed average trajectory of marital satisfaction clustered by marital status: (a) persons still in the relationship and (b) persons that ended their relationship to their married partner. The fact that the trajectories (both, for men and women) differ significantly between the marital status groups, is a strong indicator for a state-dependent TVC and thus a modelling technique able to appropriately adress this endogeneity has to be applied. Arguably rate dependence may also apply in this example, as Clark et al. (2008) found a strong anticipatory effect of the life event divorce regarding life satisfaction for men and women. It seems plausible that a similar effect exists for marital satisfaction. Besides the proper estimation in the presence of endogeneous time-varying covariates, an additional difference between the joint model and the classical time-to-event model is the inclusion of predictors for modelling the trajectory of the TVC. The basic approach of a TVC does not allow to investigate its predictors (Figure 1: Covariate 1, 2, 3 would not be taken into account for the upper panel).

Click to enlarge

Figure 2

Estimated Average Trajectory of Relationship Satisfaction of Persons by Event Status and Sex

Note. Non-linear smoother by sex (dark: men, light: women). Based on the illustration by Crowther et al. (2013).

Even though a two-stage approach can be applied to first model the TVC as a function of covariates in order to include their effects and further overcome the problem of LVCF, it has some unfavourable statistical properties mainly arising from the independent estimation of the two submodels. The subsequent time-to-event model treats the prediction of the firstly fitted longitudinal outcome, as if it was estimated without any uncertainty. This results in an underestimation of standard errors for the TVC and thus cannot be interpreted appropriately in terms of statistical inference. For an in-depth investigation on the consequences of using the aforementioned approaches instead of a joint model see Sweeting and Thompson (2011).

Compared to these two approaches (TVC approach, two-stage approach) the joint model for longitudinal and time-to-event data is advantageous as it:

A main drawback of joint models is the computational effort during the fitting procedure. In contrast, two-stage models are less computationally demanding. Therefore, some research focuses on combining the advantages of the two methods — the unbiased estimates from the joint model and the fast estimation of the two-stage approach (Leiva-Yamaguchi & Alvares, 2020).

Method

Joint models (Faucett & Thomas, 1996; Rizopoulos, 2012; Wulfsohn & Tsiatis, 1997) overcome the above mentioned shortcomings as they allow for joint modelling of a repeatedly measured outcome alongside the risk of having an event of interest. Rather than using the observations of the TVC, the joint model considers the repeatedly measured TVC as the result of a longitudinal process subject to its own model (marital satisfaction, Figure 1: upper part). This longitudinal process is combined with the related time-to-event process (risk of marriage dissolution, Figure 1: lower part).

The two submodels are described for the example in Figure 1 before merging them to the joint model. For a more general notation we refer the reader to Hickey et al. (2016).

First, the time-varying covariate is modeled via an appropriate model, in general a linear mixed model (LMM)², allowing for intra-individual variance along the time axis captured by random intercepts (b_0i) and possibly random slopes b_1i (commonly used on time as a covariate). The model corresponding to Figure 1 can be written as

with ε_i(t) ∼ N (0, σ²) where i indicates the individual and t the time point of the measurement. In vectorized form the model can be rewritten as

where $x_{i_{long}}(t)\text{'}$ is a row vector with all covariate values and a leading 1 for the intercept for person i at time t and z_i(t)′ = (1 t)′ holds the covariate values of the random effects, in this case a random intercept and a random slope on time t. Using this model allows to incorporate covariates as predictors of the estimated values of the TVC y_i(t). The covariates may be time-constant or time-varying with β being the regression coefficient vector. As in a classical LMM, the random intercepts and random slopes are assumed to stem from a multivariate normal distribution b_i ∼ $𝒩$ (0, Q).

The second related process is a time-to-event model,³ which is used to model the risk of having an event over time. The general form is a proportional hazards model, which consists of a baseline hazard h₀(t) scaled by a covariate part. The corresponding equation to Figure 1 is given by

It can be rewritten in vectorized form as

The modeled hazard function h_i(t) states the instantaneous risk of person i of having an event at time t (i.e., have a marriage dissolution). This model also contains covariates (Figure 1: Covariates 3–5), which may be exogenous time-varying or time-constant and a vector of coefficients γ.

In order to take marital satisfaction as an endogeneous covariate into the model, a joint model links the estimated value of the TVC process m_i(t) to the time-to-event model by incorporating it as a predictor and thus estimates their coefficients jointly:

The coefficient α is called the association parameter. In contrast to the two-stage approach all coefficients (β, b, γ, α) are estimated simultaneously such that all uncertainty is included in the estimation procedure. Estimating the value of marital satisfaction ${\widehat{m}}_i(t)$ involves all available observations of the person, such that the hazard in a joint model at this time does implicitly also depend on the covariate history M_i(t).

By including a covariate in both submodels, e.g., x_i3(t) in Equation (5), its direct and indirect effect on the risk of having an event can be separated. Thinking of a variable which has a strong influence on the TVC and further a smaller but significant impact on the survival: By estimating a unique $\widehat{\beta }$ coefficient as well as a $\widehat{\gamma }$ coefficient and the association parameter $\widehat{\alpha }$ we decompose the total effect via: $\widehat{\alpha }\widehat{\beta }+\widehat{\gamma }$. Hereby $\widehat{\beta }$ indicates the mediated effect of the covariate via the trajectory and $\widehat{\gamma }$ represents the direct effect on the risk of having an event. This decomposition is especially helpful to understand the effect pathway of the respective covariate.

Estimation of the coefficients can be done using a Maximum Likelihood approach (Expectation-Maximization Algorithm) or Bayesian methods (Markov-Chain Monte-Carlo sampling). The different estimation strategies for joint models are presented and compared by Rappl et al. (2021). Joint models for longitudinal and time-to-event data are implemented in several R packages (Philipson et al., 2018; Rizopoulos, 2010) and are also available in STATA (Crowther, 2020; Crowther et al., 2013). Furthermore, joint models allow for individual-specific predictions as they control for individual characteristics via the random effects in the longitudinal model.

A Different Perspective on Joint Models

In the previous section, the joint model was presented with particular emphasis on the time-to-event submodel. However, there is also a body of literature focusing on the longitudinal submodel where the time-to-event submodel is used to model informative censoring from the longitudinal study (Asar et al., 2015; Hogan & Laird, 1997; Vonesh et al., 2005; Wu & Carroll, 1988). The utilisation of a joint model may be advantageous for example in situations involving group comparisons of trajectories in the presence of censoring or missing data mechanisms related to the trajectories. One approach to account for the informative censoring is to utilise an association structure based on the random effects of the longitudinal submodel to be included in the survival submodel. This strand of literature refers to joint models as shared parameter models and contributed to the strand of literature on misspecified regression models, e.g., ignored measurement error in covariates or ignored random effects (see e.g., Bound et al., 2001; Molenberghs & Verbeke, 2001).

Data Set

The analysis is based on the German pairfam data (Brüderl et al., 2023). Pairfam (“Panel Analysis of Intimate Relationships and Family Dynamics”) is a longitudinal study with 14 annual waves contributing to shed light to changes in family and relationship structures. It started in 2008 with over 12,000 respondents. Another sample of 1,489 East-German anchor persons (“DemoDiff”) is merged as a supplementary to the data base. A detailed description of the pairfam study can be found in Huinink et al. (2011). The relationship biographies of the respondents as well as the annual questionnaire about the satisfaction with the relationship can be used to build a joint model.

The final sample consists of all persons in the pairfam data, who were married in their first marriage over the course of at least three interviews. This restriction has been made due to two reasons: First, some measurements of the longitudinal variable on marriage satisfaction are needed for proper analysis. Second, using only the first marriage of a person is based on previous empirical findings that relationship stability differs between the first and following marriages and that a selection bias may be present (Jensen et al., 2016), which might also skew the results of the model. This leaves us with a final sample size of N = 3,616 persons/marriages of which 247 (≈ 7%) stated an end of this relationship during the observation period (number of events). We did not take the actual month of divorce as event time but the stated end of relationship (marriage dissolution).

Marital satisfaction is measured as the answer on an 11-point scale to the question “All in all, how satisfied are you with your relationship?” Some exemplary trajectories of marital satisfaction are depicted in Figure 3. Note that marriages which started before the first interview such as the person in Panel C in Figure 3, are not left censored for the time-to-event model since we know the start of their marriage. They just start at a different point in time with time-dependent information. In order to allow users to reproduce the analysis, a synthesized data set can be found in the web supplementary material at Potts et al. (2026).⁴

Click to enlarge

Figure 3

Example Trajectories of Relationship Satisfaction From Pairfam Participants

Note. Vertical dashed lines indicate the time of marriage dissolution.

Table 1 summarises the time-constant variables used from the pairfam data that are included in the models. Relationship duration at marriage, age at marriage and premarital cohabitation were included as time-constant covariates and years of education, personal net income, amount of household work, labor force status, children (presence of preschool child and number of children under 18 in the respondents household) and gender role attitudes are included as time-dependent exogenous variables. The variable on the division of household work is a weighted sum index resulting from five survey items. They are measured on a 5-point scale with endpoints indicating that the respondent does all the work (high values) or the partner does all household work (low values) and thus is a relative measure between the spouses. Further, gender role attitudes is also a sum index over three items. Note that all metric variables were z-score standardized for model estimation purposes and the scale of the time is changed from months to vary in the interval [0, 1] with 1 being the overall latest time-point measured in a marriage in the sample. Further note that for some of the covariates labelled as exogenous, the exogeneity assumption can also be called into question, as they may also respond to couple-separation (e.g., labor market status, number of children). However, for the sake of simplicity, we focus on one single endogenous covariate in this tutorial (for possible extensions see the Summary & Conclusion section).

Model Specification

Based on the reviewed literature on marital satisfaction and marital dissolution, two joint models will be fitted separately for men and women. They are identical in terms of included covariates, method and association structure.

The longitudinal model on satisfaction with the marriage will be modelled by an LMM including a random intercept and a random slope term for the duration of marriage (t). We included a random slope since intercepts and slopes differ between individuals (see Figure 3 for some example trajectories). The above mentioned variables are taken into account as covariates to model the longitudinal variable properly. Hereby, i is an indicator for the person and j denotes the measurement at time-point j:

The time-to-event model on time to marriage dissolution will be modelled jointly with the longitudinal model from Equation (6), where the estimated values of satisfaction explain the risk of marriage dissolution. Specifically, the time-to-event model includes a B-Spline approximation of the baseline hazard h₀(t).

Equation (7) shows the final representation of our main joint model on marital satisfaction and the risk of marriage dissolution.

Implementation

For implementation of the above model we use the package $\mathrm{JM}$ (Rizopoulos, 2010) for $R$ (Version 4.4.0) (R Core Team, 2024). $\mathrm{JM}$ combines two models that are built with their specific $R$ packages. The longitudinal model is constructed using $\mathrm{lme()}$ from the $\mathrm{nlme}$ package (Pinheiro et al., 2023) and therefore requires the usual data structure in long format, where each individual spans several rows corresponding to the observation time points, each holding the covariate value of the time point, respectively.

The time-to-event model is fitted using $\mathrm{coxph()}$ from $\mathrm{survival}$ (Therneau, 2024). The structure of the underlying data set is equivalent to the long format start-stop-event logic when further exogenous TVCs are used, i.e., several rows per individual, indicating the current measured values. In case of lack of other exogenous TVCs in the model, the data set for the time-to-event model reduces to one row per individual and supplements the long format data set for the longitudinal submodel (see Rizopoulos (2012)). Since our model includes other exogenous TVCs (e.g., labor force status) one single data set in the classical start-stop-event logic is used.

Both models serve as inputs for the final $\mathrm{jointModel()}$ command. For further information on the (optional) arguments in the $\mathrm{jointModel()}$ function, we refer to Rizopoulos (2010). $R$-code and a synthesized data set to replicate the example can be found in the web appendix of this paper (see Potts et al., 2026).

# time-to-event model
	modsurv_female <- coxph(Surv(time = t1, time2 = t2, event = status)~
	  yeduc + ageatm + preschoolchild + nchild + premarcohab +
	  sw_weight + incnet + relduratmar + lfs_rec + genderatt_s,
	  data = df_female, x = TRUE, model=T, cluster = id)

# longitudinal model
	modlong_female <- lme(sat31 ~ t + I(t^2) + sw_weight + ageatm +
	  preschoolchild + nchild + premarcohab + yeduc + incnet +
	  relduratmar + lfs_rec + genderatt_s,
	  data = df_female, random = ~ t | id)

# joint model for longitudinal and time-to-event data
	modjoint_female <- jointModel(modlong_female, modsurv_female,
	  timeVar = "t", method = "spline-PH-GH",
	  control = list(verbose=T, iter.EM=100))

Estimation Results

Separate models were fitted for men and women. This section starts with the estimation results of the joint model for women.

The longitudinal submodel estimates the relationship between the covariates and marital satisfaction (outcome) (Table 2, left side). The model results in a U-shaped effect of time. The number of years of education shows a positive, significant effect on marital satisfaction. Some other covariates show negative, linear effects on satisfaction with the relationship: higher age at marriage and higher personal net income are associated with a significantly lower marital satisfaction. Women with children (compared to childless women) also reveal lower values of satisfaction. There is even an additional negative effect if there are preschool children present in the household. The index for the division of household work reveals a negative and statistically significant effect, i.e., women who stated to do more household work are less satisfied with their relationship during marriage. The relationship duration at time of marriage and the gender role attitudes show insignificant coefficients in this model. Next, we examine the time-to-event submodel for women (Table 2, right side) including the properly modelled endogenous variable marital satisfaction. This model indicates which variables still have a direct effect on the risk of marital dissolution when controlling for marital satisfaction. Looking at the association parameter (last row in Table 2), we observe the expected strong negative effect of marital satisfaction on the risk of marriage dissolution: the higher the current value of satisfaction with the relationship, the lower the risk of marriage dissolution. Besides this relationship, there are only few significant effects in the time-to-event submodel. Higher educated women and women that were in a long relationship with their married partner before marriage have a lower risk of marriage dissolution. As an example for the decomposition of effects, we focus on the variable of relative household work in the following. In contrast to other research results (e.g., Frisco & Williams, 2003)⁵, the relative load of household work done by a person has no direct effect for women on their risk of marriage dissolution. There still remains the indirect effect via the mediator of marital satisfaction: a higher share of household work done by the female respondent results in a significantly lower marital satisfaction which results in significantly higher risk of marriage dissolution. The estimated total effect of the index on the divison of household work on the risk of marriage dissolution adds up to −0.5552 × −0.1369 + 0.1089 ≈ 0.1849.

For men, the results reveal interesting differences in both submodels (Table 3): Regarding the marital satisfaction (longitudinal model ), socio-economic variables such as education as well as personal net income do not show significant effects. Similarly to the model for women, children decrease the relationship satisfaction in the marriage compared to childless men and the effect size is larger than for women. In contrast to the model for women, there is no additional significant, negative effect of preschool children. Even though the labor force status shows an influence on the marital satisfaction for men, these results have to be interpreted with caution, as over 80% of the person periods for men indicate a full-time employment. The relative load of household work has a smaller effect on marital satisfaction and is also negative and statistically significant.

The estimated association parameter $\widehat{\alpha }$ in the time-to-event submodel is also significant but the effect size is smaller compared to the model for women. In other words, the estimated marital satisfaction appears less predictive for the risk of marriage dissolution for men than for women. Another difference between the sexes in the submodel regarding the risk of marriage dissolution is the effect of the division of household work: This covariate has a positive and statistically significant effect on the risk of marriage dissolution in the model on men (p = 0.0660). In contrast to the model for women, the indirect mediated effect is supplemented by a direct effect of household work on the risk of marriage dissolution. The estimated total effect of this predictor variable on the risk of marriage dissolution can be derived as −0.4534 × −0.0603 + 0.2413 ≈ 0.2686.

Besides the decomposition of interpretable covariate effects, another strength of a joint model is the option to predict individual probabilities to still be in the relationship after the last measurement of a person. This may be useful for intervention planning in different research questions on the micro-level. In contrast to a classical TVC approach, these predictions are based on the whole estimated longitudinal trajectory of marriage satisfaction and do not only rely on the last measured value. The orange line in Figure 4 shows the predictions of a fictional part-time working woman, who lived together with her partner before marriage. She lives with two children, at least one being a preschool child and shows median values (for women) for the other covariates. The only part that varies between the four plots is her trajectory of marital satisfaction indicated by the crosses left to the dashed line showing the time of the last measurement. Looking at a certain point of time in the future (black solid line), the different trajectories result in different predicted probabilities to still be in the marriage, even though each trajectory ends with the same last measurement of a 4 on the

Click to enlarge

Figure 4

Predicted Probability of Still Being in the Marriage for a Person Varying Only the Marital Satisfaction Trajectory

Note. Covariate values: part-time working, 2 children, at least one preschool child, premarital cohabitation, median values (for women) for the other covariates.

11-point scale. This feature also allows to update the prediction with every newly obtained measurement on a person such that a researcher can trace the development dynamically (see Appendix A.4).

Comparison to Other Modelling Approaches

In this section, we compare the estimation results, i.e., coefficients, standard errors, predictive performance (Mean Squared Error), of the joint model with two other modelling approaches (TVC approach in a Cox-proportional hazards model and a two-stage model). Comparison tables of the estimation results can be found in Appendix A.1.

A classical TVC approach in a time-to-event model ⁶ would heavily underestimate the effect of marital satisfaction on the risk of marriage dissolution (e.g., women: TVC approach: −0.3083, joint model: −0.5552). This may result from the fact that the risk of marriage dissolution is not only dependent on the observed current value but due to the estimation routine on the whole trajectory of the marriage satisfaction until the point in time. Modelling marital satisfaction as an endogenous covariate results in some differences regarding the other covariates in the time-to-event submodel as well. For example, the effect of the division of household work for women on the risk of marriage dissolution is overestimated in a classical TVC model and the standard error is underestimated which leads to a smaller p-value. These differences highlight the importance of the correct model choice when dealing with endogenous covariates in time-to-event models.

In our example the differences between a two-stage model and a joint model are only small and the tendency to underestimate the uncertainty of the satisfaction variable can only be found in the model for women. The regression coefficients for marital satisfaction differ only slightly with −0.5347 in the two-stage model and −0.5552 in the joint model, and the standard error of 0.0530 in the two-stage model is smaller than the standard error of 0.0551 of the presented joint model (men: two-stage: −0.4240, s.e. = 0.0708; joint model: −0.4534, s.e. = 0.0702).

Focussing on the longitudinal submodel, we compare the standalone longitudinal model (mixed model fitted with $\mathrm{lme())}$ and the longitudinal submodel of a joint model which controls for the non-random drop-out due to marriage dissolution. There are only small differences in coefficients and standard errors between the two models. However, in the joint model for women the effect of time of marriage has a more pronounced U-shape i.e., larger absolute coefficients for the linear and quadratic term. In the model for men, modelling the non-random drop-out by marriage dissolution leads to a change in the size and significance of the coefficient of premarital cohabitation. Note that the differences between the longitudinal model approaches may be larger in other applications when a larger share of events (higher number of drop-outs, i.e., more marriage dissolutions) is present.

We further evaluated the predictive performance of the three modelling approaches via predicting the event probability for several time points after the last individual measurement of marital satisfaction. Figure A.1 shows the mean squared error (MSE) of the models at specific points in time and highlights the advantages of a joint model regarding predictions in the longrun. While the MSE is smallest for the TVC approach when predicting the event outcome up until six months after the last measurement, the MSE for the joint model is smaller when looking at later times and outperforms the two competing models after ten months.