Flexible parametric multistate modelling of employment history

Abstract

A multistate model is used to describe employment history. Transition-specific rates are defined using generalized gamma distributions and Gompertz distributions. This flexible parametric modelling of the rate of change is combined with latent classes for unobserved propensity to change jobs. The propensity is described by two latent classes which can be interpreted as consisting of movers and stayers. The modelling is illustrated by analysing longitudinal data from the German Life History Study.

Keywords

demography frailty model generalized gamma distribution Gompertz distribution Markov process survival analysis

1 Introduction

Multistate models are used to describe stochastic processes where the change of status is of interest. Many applications can be found in biostatistics with the illness-death model as the quintessential example. An illness-death model can be defined by three states: a healthy state, an ill state and the dead state. Of interest can be risk factors for the onset of the illness or expected duration in the ill state. In social statistics and in demography, multistate models are used to study processes such as changes in region of residence, employment history, or changes in marital status.

Statistical methods for multistate models are typically not discipline-specific. As an example, a three-state illness-death model for medical data can be quite similar to a three-state model for employment history when the states in the latter are defined as employed, unemployed and retired.

The aim of this article is to introduce the generalized gamma distribution (Stacy, 1962) for flexible parametric continuous-time multistate modelling in demography. The application is with respect to employment history, where states are defined corresponding to the number of past jobs and intermediate spells without a job. The time scale in our model for employment history is time since entry into the labour market. We use the generalized gamma distribution to model transition-specific time dependence, and we show that using this distribution improves the statistical inference when combined with transition-specific exponential and Gompertz distributions. In addition, we show that parametric models can be extended with a definition of latent classes.

For continuous-time multistate models, the generalized gamma distribution is discussed in Jackson (2016). The parametric modelling is within a wider methodological framework for multistate models described by, for example, Kalbfleisch and Lawless (1985), Kay (1986), Hougaard (2000), Jackson (2011), and Van den Hout (2017). The latent-class model that we define can be seen as a random-effects model with a discrete distribution for the random effects. Latent-class models for discrete-time multistate models are discussed in, for example, Vermunt et al. (1999) and Bartolucci et al. (2012). An example in demography is Dias and Willekens (2005), who discussed determining the number of latent classes. For continuous-time multistate models which include random effects see, for example, Hougaard (2000), Putter and Van Houwelingen (2015), and Van den Hout (2017). For multistate modelling in demography, an overview with a wide range of random-effects structures is given by Bijwaard (2014).

When multistate models are applied in biostatistics, we can see the latent-class model as a frailty model. For example, if progression through a set of states denotes a deterioration of health, a model with two latent classes can distinguish individuals who move quickly through the states (the frail individuals) from those who move less quickly (the more healthy ones). In the current article, the latent classes are defined with respect to employment history and will allow us to distinguish individuals who tend to change job more quickly (the movers) from those who tend to stay put (the stayers). In the application, we will show that such a distinction can lead to a model that fits the data better.

Our multistate model for employment history combines the generalized gamma distribution with the latent-class approach. Although this distribution and the latent-class approach have been discussed in the literature separately (see the references above), our contribution is to combine these two concepts and thus defining a very flexible statistical modelling framework. The article illustrates this by defining a series of models, and by model comparison, validation and interpretation. We show that the modelling framework is general and allows for a wide range of applications in demography.

2. German Life History Survey

The German Life History Study (1980–2005) provides retrospective life course information for Germans born between 1919 and 1951 (Mayer, 2015). The study is often used to investigate education, employment history and family formation.

For the current article, we use the Blossfeld–Rohwer subsample of the German Life History Study (Blossfeld and Rohwer, 2002). These data are available in the Biograph package in the R software. This package is introduced by Willekens (2014), who used the Blossfeld–Rohwer subsample in a three-state model for employment history.

Figure 1

Eleven-state process for employment history in GLHS

This subsample contains data for 201 individuals on job episodes and spells without a job. The start and end of each job episode are available on a time scale in months. Individual background information is given by covariates such as age, gender and education.

Interviews for GLHS were conducted in 1981. In what follows, we propose a statistical model for employment history as known at the time of the interview. This history is defined by a series of mutually exclusive states. Because not many individuals in GLHS have had more than six jobs, we restrict the modelling of employment history up to (and including) the sixth job.

Figure 1 shows the process that we will model in this article. States 1–6 denote the first job up to the sixth job. States 7–11 denote episodes between jobs; that is, state 7 denotes having no job after the first job, state 8 denotes having no job after the second job, and so on. With the imposed restriction to a maximum of six jobs, state 6 acts as an absorbing state.

It is possible to define a process with fewer states by collapsing the five no-job states into one state. However, the model comparison in Section 6 shows that it is worthwhile to distinguish the five no-job states, as illustrated in Figure 1.

In the subsample (denoted GLHS from now on), all individuals start in state 1, but they do so at different times and different ages. The time scale in GLHS is century months. For example, the first individual in the data is a man who starts his first job at century month 555. In calendar years, this is $(555 - 1) / 12 + 1900 = 1946$ . In GLHS, three birth cohorts are represented: 1929–1931, 1939–1941, and 1949–1951.

For the 201 individuals in GLHS, Figure 2 shows the start time of the first job (in calendar years) and the age at which the individuals started their first job. The graph shows a cohort effect in the sense that the age at the start of the first job is increasing with increasing calendar year for the start of the first job. Therefore, we will include cohort effects in the statistical model for employment history.

Table 1 is the state table for the 11-state employment history in GLHS. This table shows the frequencies for successive pair of states. Right-censored histories are presented in the $(r, r)$ entries with $r = 1, . . ., 5,7, . . ., 11$ . To give an example, of the 201 individuals who started in state 1, 16 are still in their first jobs at the time of the interview, 96 went straight from their first jobs to their second job, and 89 had a period without a job after their first job. The table shows that data are sparse for states 5, 6, 10 and 11. In totals there, are 564 job episodes (sum of the frequencies in the first five rows) and 159 spells without a job (sum of the frequencies in the last five rows).

Figure 2

For the 201 individuals in GLHS, a two-dimensional representation of the start time of the first job and the age at the start of the first job

Table 1

State table for GLHS: number of times for each successive pair of states. The states are defined by employment history; see Figure 1

	To
From	1	2	3	4	5	6	7	8	9	10	11
1	16	96	0	0	0	0	89	0	0	0	0
2	0	44	78	0	0	0	0	40	0	0	0
3	0	0	39	52	0	0	0	0	16	0	0
4	0	0	0	27	24	0	0	0	0	11	0
5	0	0	0	0	11	18	0	0	0	0	3
7	0	66	0	0	0	0	23	0	0	0	0
8	0	0	29	0	0	0	0	11	0	0	0
9	0	0	0	10	0	0	0	0	6	0	0
10	0	0	0	0	8	0	0	0	0	3	0
11	0	0	0	0	0	2	0	0	0	0	1

To model potential time dependence in employment history, we need to specify a time scale for the process. Possible options are time in current job, age and time since entry into the labour market. For our statistical model, we use the third option following the choice of Blossfeld and Rohwer (2002, Chapter2) for their model for the rate of leaving the current job. This option implies, for example, that our modelling allows two individuals who have been in their third job equally long to have a different distribution for moving to the fourth job due to a difference in time since entry into the labour market. It allows us to investigate the effects of gender and education while controlling for time spent in the labour market.

3 Distributions for transition times

To model change of state in continuous time, we specify continuous parametric distributions for transition time $T \geq 0$ . This section starts with the definition of transition-specific hazard functions and state-specific survivor functions and will then discuss the exponential distribution, the Gompertz distribution and the generalized gamma distribution.

Our modelling is based on continuous-time stochastic processes, which are defined as ${Y_{t} | t \in (0, \infty)}$ with $Y_{t}$ denoting the state at time t. The range of $Y_{t}$ is defined by the state space $S$ .

In our application, we have finite $S = {1,2, . . ., 11}$ , with state 6 being the only absorbing state. Let $(t_{1}, t_{2}]$ denote a generic time interval, and define $Y_{j} = Y_{t_{j}}$ , for $j = 1,2$ . Specific for GLHS, the following type of observations are possible for a time interval $(t_{1}, t_{2}]$ . If there is no transition, then $(y_{1}, y_{2}) = (r, r)$ , for $r \in a l S ∖ {6}$ . If there is a transition, then $(y_{1}, y_{2}) = (r, s)$ with $s = r + 1$ or $s = r + 6$ if $r \in {1,2, 3,4, 5}$ , and $s = r - 5$ for $r \in {7,8, 9,10,11}$ . For example, if the individual moves out of state 1, then he or she moves to state 2 or state 7. Given state 7, the next possible state is 2.

Let $h_{r s} (t)$ denote the hazard for moving from state r to state s at time t; that is,

h_{r s} (t) = \lim_{Δ ↓ 0} \frac{P (Y_{t + Δ} = s | Y_{t} = r)}{Δ} .

If a transition is not possible from r to s, we define $h_{r s} (t) = 0$ . Let $h_{r •} (t)$ denote the hazard for moving out of state r. For example, for state 1 we have $h_{1 •} (t) = h_{12} (t) + h_{17} (t)$ .

For $(t_{1}, t_{2}]$ , the conditional survivor function for state r at $t_{1}$ is given by

S_{r} (t_{2} | t_{1}) = \exp (- \int_{t_{1}}^{t_{2}} h_{r •} (u) d u) .

This definition of $S_{r} (t_{2} | t_{1})$ implies that given state r at $t_{1}$ , the hazard function $h_{r •}$ defines the probability of a transition out of state r. In the application, if all five hazard functions are defined by exponential distributions, then the corresponding transition probabilities define a time-homogeneous first-order Markov chain. If the hazard functions are time-dependent or contain covariate effects, then the process is not a Markov chain in the strict sense; see also Commenges (1999).

3.1 Exponential and Gompertz distribution

The transition-specific hazard function for the exponential distribution is given $h_{r s} (t) = \exp (μ_{r s})$ , where $μ_{r s} \in R$ . The exponential model is a very common model for multistate processes and is readily available in software for multistate models. Note that this hazard function does not depend on time t, which may be too restrictive in some applications.

The transition-specific hazard function for the Gompertz distribution is $h_{r} (t) = \exp (μ_{r s} + ξ_{r s} t)$ , where $μ_{r s}, ξ_{r s} \in R$ . As a distribution for a time to event, the Gompertz distribution is used in biostatistics and sociology; see, for example, Hougaard (2000) and Blossfeld and Rohwer (2002), respectively. For usage in demography; see, for example, Mueller et al. (1995), Olshansky and Carnes (1997), and Bongaarts (2005)

Models with exponential and Gompertz distributions can be extended in the usual way by log-linear regression: $μ_{r s} = β_{r s}^{⊤} x$ , where $x$ is a vector of covariates values including an intercept.

3.2 Generalized gamma distribution

The generalized gamma distribution is a very flexible continuous parametric distribution for event time $T \geq 0$ . The distribution includes the exponential, Weibull, log-normal and gamma distributions as special cases; that is, these cases can be definedby restricting certain parameters in the generalized gamma distribution. We will first review the generalized gamma distribution and then use this distribution to define the transition-specific hazards.

Following the presentation and terminology in Cox et al. (2007), the generalized gamma distribution $T \sim G G (μ, σ, q)$ has probability density function

f (t | μ, σ, q) = \frac{| q |}{σ t Γ (q^{- 2})} {(q^{- 2} (\exp (- μ) t)^{q / σ})}^{q^{- 2}} \exp (- q^{- 2} (\exp (- μ) t)^{q / σ}),

3.1

for location μ, scale $σ > 0$ , and shape q. The survivor function is given by

S (t) = {\begin{array}{l} 1 - Γ [q^{- 2} (\exp (- μ) t)^{q / σ}; q^{- 2}] & i f & q > 0 \\ Γ [q^{- 2} (\exp (- μ) t)^{q / σ}; q^{- 2}] & i f & q < 0, \end{array}

3.2

where $Γ [t; γ] = \int_{0}^{t} x^{γ - 1} e^{- x} d x / Γ (γ)$ . Prentice (1974) introduced the above parameterization for the generalized gamma by Stacy (1962). This choice of parameterization allows $q \leq 0$ . The generalized gamma distribution is available in R via the flexsurv package (Jackson, 2016).

The log-normal distribution is defined as $G G (μ, σ, q = 0)$ ; the Weibull distribution as $G G (μ, σ, q = 1)$ and the gamma distribution as $G G (μ, σ, q = σ)$ . To illustrate this for the Weibull distribution, for $q = 1$ density (3.1) reduces to

f (t | μ, σ, q = 1) = \frac{1}{σ} t^{1 / σ - 1} (\exp (- μ))^{1 / σ} \exp (- (\exp (- μ) t)^{1 / σ}) .

Re-parameterization by $τ = 1 / σ$ and $λ = \exp (μ)$ gives

f (t | τ, λ) = \frac{τ}{λ} {(\frac{t}{λ})}^{τ - 1} \exp (- {(\frac{t}{λ})}^{τ}),

which is a standard representation of the Weibull density. The median of this distribution is $λ {(\log (2))}^{1 / τ} = \exp (μ) {(\log (2))}^{σ}$ . This illustrates that the location parameter μ determines the median, given a fixed value for scale $σ$ . This also holds for the unrestricted generalized gamma distribution. Scale $σ$ and shape q together determine the type of hazard function as stated and illustrated in Cox et al. (2007); see also the data analysis in Section 6.

For the multistate model, we define transition-specific hazard as $h_{r s} (t) = f_{r s} (t) / S_{r s} (t)$ , where the density function $f_{r s} (t)$ and the survivor function $S_{r s} (t)$ are defined by the generalized gamma distribution $G G (μ_{r s}, σ_{r s}, q_{r s})$ . If there are two possible transitions out of state r, then we are in a competing-risks setting. For example, the function $S_{12} (t)$ as defined above should not be interpreted as the survivor function with respect to being in state 1; see also (Putter et al., 2007, Section 3.2).

4 Latent-class models

A latent-class model is defined by specifying K classes, probabilities $π_{i k} = P (i \in c l a s s k)$ for individual $i \in {1, . . ., N}$ and class $k \in {1,2, . . ., K}$ and class-specific model parameters. The probabilities $π_{i k}$ model the class membership with the restriction $\sum_{k = 1}^{K} π_{i k} = 1$ . Class-specific model parameters can refer to location, shape or scale parameters that vary across the K classes.

The class-specific parameters induce class-specific distributions for the transition times. As an example, consider the exponential model and $K = 2$ . For instance, we want to explore the latent-class modelling for direct transitions from one job to the next job. We can define $π_{i 1} = π_{1}$ and $π_{i 2} = 1 - ð_{1}$ , for all i, and make the relevant location parameters class-specific by defining $μ_{r s . k} = β_{r s} + b_{k}$ , where $(r, s) \in {(1,2), (2,3), (3,4), (4,5), (5,6)}$ and $b_{1} = - b_{2}$ . The two latent classes induce two extra parameters to be estimated: $π_{1}$ and $b_{1}$ . The first is the probability to belong to class 1, and the second is the effect on the location parameter due to being in class 1. For example, if $b_{1} < 0$ , then $\exp (μ_{r s . 1}) < \exp (μ_{r s . 2})$ , which implies that individuals in class 2 are more likely to move to state s than individuals in class 1. A model like this allows us to distinguish so-called movers from stayers: unobserved heterogeneity in the propensity to move to a next state.

The distinction between individuals who move quickly through the states and individuals who tend not to move so quickly is described as tracking by Satten (1999): there is a correlation between the transition times within individuals, and rapid progressors can be distinguished from slow progressors. In the literature, a mover/stayer model can also denote a stochastic process that is a mixture of two processes, one of which has transition hazards equal to zero; see, for example, Frydman (1984) and Cole et al. (2005). This option, however, will not be explored in the current article.

Latent-class models can be seen as random-effects models where the random effects have a discrete distribution; see, for example, the generalized linear mixed models as discussed in Aitkin (1999). When fitting the latent-class as a random-effects model, the classes that are created are primarily aimed at capturing unobserved heterogeneity interpretation of the classes is not of primary interest.

5 Estimation

For individual i, for $i \in {1,2, . . ., N}$ , GLHS data were formatted such that the data are given by times $t_{i} = (t_{i 1}, t_{i 2}, . . ., t_{i n_{i}})$ , states $y_{i} = (y_{i 1}, y_{i 2}, . . ., y_{i n_{i}})$ and baseline covariates values $x_{i} = (x_{i 1}, x_{i 2}, . . ., x_{i m}$ ). We also define transition indicators $δ_{i} = (δ_{i 2}, δ_{i 3}, . . ., δ_{i n_{i}})$ to maintain the standard notation in survival analysis. Note, however, that in the GLHS data, times $t_{i 1}, . . ., t_{i n_{i} - 1}$ are all transition times, and $t_{i n_{i}}$ is either a right-censored time, or a transition time for moving to state 6.

Consider a time interval $(t_{i j}, t_{i j + 1}]$ , $j = 1, . . ., n_{i} - 1$ . When $y_{i j} = y_{i j + 1}$ , there is no transition, and the conditional survivor function for observing state $Y_{i j + 1}$ at $t_{i j + 1}$ is $S_{y_{i j}} (t_{i j + 1} | t_{i j})$ . When $y_{i j} \neq y_{i j + 1}$ , time $t_{i j + 1}$ denotes a time of a transition, and the density function for moving to observed state $Y_{i j + 1}$ at $t_{i j + 1}$ is $S_{y_{i j}} (t_{i j + 1} | t_{i j}) h_{y_{i j} y_{i j + 1}} (t_{i j + 1})$ . The former case is labelled by $t a_{i j + 1} = 0$ and the latter by $δ_{i j + 1} = 1$ .

Individual contributions to the likelihood function for the fixed-effects models are given by

L_{i} (θ | t_{i}, δ_{i}, y_{i}, x_{i}) = \prod_{j = 1}^{n_{i} - 1} S_{y_{i j}} (t_{i j + 1} | t_{i j}) h_{y_{i j} y_{i j + 1}} {(t_{i j + 1})}^{(δ_{i j + 1})},

for $i = 1,2, . . ., N$ , with $θ$ the vector containing all the model parameters. The overall likelihood function is $L (θ) = \prod_{i = 1}^{N} L_{i} (θ | t_{i}, δ_{i}, y_{i}, x_{i})$ .

Individual contributions to the likelihood function for the latent-class models are given by

L_{i} (θ | t_{i}, δ_{i}, y_{i}, x_{i}) = \sum_{k = 1}^{K} (\prod_{j = 1}^{n_{i} - 1} S_{y_{i j}} (t_{i j + 1} | t_{i j}, k) h_{y_{i j} y_{i j + 1}} {(t_{i j + 1} | t_{i j}, k)}^{(δ_{i j + 1})}) π_{i k},

where $θ$ includes the parameters for the latent-class modelling. Probabilities $π_{i k}$ can be linked to individual covariate values using a multinomial regression model. However, in the GLHS data analysis in Section 1, weassume $π_{i k} = π_{k} \forall i$ .

The likelihood function for the fixed-effects models is the standard format for survival data extended to multiple event times within individuals. Note that possible association between observations within an individual is not modelled. In contrast, the likelihood function for the latent-class models allows a correlation between observations through the latent-class parameters.

For the maximum likelihood estimation, we use the general-purpose optimizer in the R software optim (R Core Team, 2013). This versatile optimizer can maximize the log-likelihood function without the need to provide explicit expressions for the derivatives. The optimizer will return the numerically differentiated Hessian matrix if requested.

Covariance of a function of estimated model parameters can be derived by using simulation. An important example of such a function are transition probabilities for a specified time interval. For the simulation, parameter vectors $θ^{(s)}$ , for $s = 1, . . ., S$ , are drawn from $N (\hat{θ}, {\hat{V}}_{θ})$ , where ${\hat{V}}_{θ}$ is the estimated covariance matrix derived from the estimated Hessian. Next, for any function $g (θ)$ , summary statistics such as mean and covariance are derived from $g (θ^{(1)}), . . ., g (θ^{(S)})$ ; see Mandel (2013) for a justification of this method.

6 Data analysis

We start the data analysis for GLHS with the basic exponential model. This model can be seen as an intercepts-only model for the location parameters and has as many parameters as transitions in the process: $μ_{r s} = β_{r s . 0}$ , for the pairs $(r, s)$ as defined by the process in Figure 1. This model has Akaike information criterion (AIC) equal to 6605.4.

Dummy variables for the cohort effects are defined as $(x_{1}, x_{2}) = (0,0)$ for birth cohort 1929–1931, $(x_{1}, x_{2}) = (1,0)$ for cohort 1939–41, and $(x_{1}, x_{2}) = (0,1)$ for cohort 1949–51; see also Blossfeld and Rohwer (2002). Let $x_{3}$ and $x_{4}$ denote the dummy variables for the effect of gender and education. Women are coded by $x_{3} = 0$ , men by $x_{3} = 1$ . Education is defined by years of education before entry into labour market: $x_{4} = 0$ for less than 12 years, $x_{4} = 1$ for 12 years or more. The threshold 12 years is equal to middle school qualification with vocational training.

For the modelling of the effects of the dummy variables, we distinguish three types of transitions: from one job to the next (A), from a job to no job (B) and from no job to the next job (C). The model for the location parameters with the covariates is thus given by

\begin{array}{l} A) μ_{r s} = β_{r s .0} + β_{A .1} x_{1} + β_{A .2} x_{2} + β_{A .3} x_{3} + β_{A .4} x_{4} {\begin{array}{l} From a job to the next : \\ (r, s) = (r, r + 1), r \in {1, ..., 5} \end{array} \\ B) μ_{r s} = β_{r s .0} + β_{B .1} x_{1} + β_{B .2} x_{2} + β_{B .3} x_{3} + β_{B .4} x_{4} {\begin{array}{l} From a job to no job : \\ (r, s) = (r, r + 6), r \in {1, ..., 5} \end{array} \\ C) μ_{r s} = β_{r s .0} + β_{C .1} x_{1} + β_{C .2} x_{2} + β_{C .3} x_{3} + β_{C .4} x_{4} {\begin{array}{l} From no job to the next job : \\ (r, s) = (r, r - 5), r \in {7, ..., 11} \end{array} \end{array}

6.1

and has AIC = 6437.2, which shows a clear improvement over the intercepts-only model.

It is possible to restrict the intercepts $β_{r s . 0}$ in such a way that the hazards for moving to a no-job state are the same: $β_{r, r + 6.0} = β_{B . 0}$ for $r \in {1, . . ., 5}$ , and the hazards for leaving a no-job state are the same: $β_{r, r - 5.0} = β_{C . 0}$ for $r \in {7, . . ., 11}$ . This defines 7 independent intercepts. Adding the 12 parameters for the dummy variables, the resulting model has AIC = 6453.9. The increase of the AIC compared to the previous model shows that it is worthwhile to distinguish the transition hazards for the five no-job states. This implies that we should not merge the five no-job states into one general no-job state.

To investigate non-constant transition hazards, we discuss two models. For the first model, we add Gompertz $ξ$ -parameters to the first four transitions from one job to the next; that is, we add $ξ_{r s}$ for $(r, s) = (r, r + 1)$ , $r \in {1,2, 3,4}$ . Because data are scarce for the transition from state 5 to state 6, we do not extend the model for this transition. The Gompertz model yields an AIC = 6400.9, which shows that allowing hazards to change over time can lead to a model with better fit.

Table 1 shows that there is good information for moving out of states 1 and 2. To allow for more flexible parametric shapes, we fit generalized gamma distributions for $(r, s) \in {(1,2), (1,7)}$ . The shape parameters are $q_{12}$ and $q_{17}$ , and the scale parameters are

σ_{12} = \exp (ψ_{12}) and σ_{17} = \exp (ψ_{17}) .

For $(r, s) \in {(2,3), (2,8), (3,4), (4,5)}$ , we use Gompertz distributions. Note that this also implies a non-constant hazard for moving from state 2 to no-job state 8. For the remaining transitions, we fit exponential distributions. The covariate effects are again as specified in (6.1). We call this Model $A$ . It has AIC = 6379.6 and shows a further improvement over the earlier models.

Next, we extend Model $A$ with two latent classes ( $K = 2$ ). We define class-specific distributions via the location parameters for the first four direct job changes. To impose the mover/stayer distinction, the class effect for the four job changes is assumed to be the same. For $k \in {1,2}$ , this implies the following change for the location parameters:

μ_{r s . k} = β_{r s .0} + b_{k} + β_{A .1} x_{1} + β_{A .2} x_{2} + β_{A .3} x_{3} + β_{A .4} x_{4},

6.2

for $(r, s) = (r, r + 1), r \in {1, . . ., 4}$ . The model is identified by the sum-to-zero restriction for the class-specific effects; that is, $b_{1} + b_{2} = 0$ . We call this Model $B$ . It has two more parameters ( $π_{1}$ and $b_{1}$ ) compared to Model $A$ and has AIC = 6366.8.

Model $B$ is a better model than $A$ , given the improvement in AIC. In the remainder of this section, we present statistical inference for Model $B$ .

Estimated parameters are presented in Table 2. The probability

π_{1}

to be in class 1 is estimated at 0.276 (with estimated standard error 0.083). For this class, the additive effect

b_{1}

on the location parameters in (6.2) is 0.641 (0.108). Given that

b_{1} > 0

, class 1 represents the movers those individuals who are more likely to move directly to a next job when compared to individuals in class 2.

Table 2

Parameter estimates for the mover/stayer model for the GLHS data on employment history (Model $B$ ). Time scale is months since entry into the labour market. Estimated standard errors in parentheses. Estimated intercepts are not included

Location				Scale and shape
$β_{A . 1}$	0.417 (0.201)	$β_{A . 3}$	0.200 (0.160)	$ψ_{12}$	0.044 (0.197)
$β_{B . 1}$	0.497 (0.205)	$β_{B . 3}$	$- 0.631$ (0.177)	$ψ_{17}$	0.311 (0.086)
$β_{C . 1}$	0.539 (0.231)	$β_{C . 3}$	2.107 (0.202)	$q_{12}$	$-$ 2.038 (0.760)
$β_{A . 2}$	0.228 (0.211)	$β_{A . 4}$	0.003 (0.188)	$q_{17}$	$-$ 0.439 (0.459)
$β_{B . 2}$	0.599 (0.206)	$β_{B . 4}$	$- 0.160$ (0.202)
$β_{C . 2}$	0.429 (0.254)	$β_{C . 4}$	$- 0.717$ (0.262)	$ξ_{23}$	$-$ 0.004 (0.002)
				$ξ_{28}$	$-$ 0.007 (0.002)
Class 1				$ξ_{34}$	$-$ 0.003 (0.002)
$π_{1}$	0.276 (0.083)	$b_{1}$	0.641 (0.108)	$ξ_{45}$	$-$ 0.004 (0.002)

For women born in 1939–1941 and with less than 12 years of education, Figure 3 depicts the hazards in Model $B$ that are time-dependent. The difference between the two classes is illustrated with the black lines representing class 1, and the grey representing class 2. Confidence bands in this figure are derived by simulation; see end of Section 5. For $1 \to 7$ and $2 \to 8$ , the hazards are restricted to be the same for both classes. Hazards are modelled conditional on current state. For example, if the women are in state 3 quickly after having entered the labour market, then the model implies a high hazard for moving the state 4.

The flexibility of the generalized gamma distribution is illustrated in Figure 3 by the arc-shape hazards transitions $1 \to 2$ and $1 \to 7$ . For example, the hazard for moving directly to the second job increases in the first 4–5 years and shows a decrease afterwards. This can be linked to the estimated scale ${\hat{σ}}_{12} = \exp ({\hat{ψ}}_{12})$ and shape ${\hat{q}}_{12} < 0$ in Table 2. This combination of values is in the set ${(σ, q) | q < \min {σ, 1 / σ}}$ that defines arc-shape hazards, starting and ending at zero Cox et al. (2007). For the Gompertz model, the estimated $ξ$ -parameters for transitions $2 \to 3$ , $2 \to 8$ , $3 \to 4$ and $4 \to 5$ are negative and illustrate a negative time effect on the hazard.

Looking at the estimated standard errors in Table 1, we see that there are clear cohort effects. The younger cohorts have higher hazards for all transitions. The effects for gender illustrate that men (coded by $x_{3} = 1$ ) are more likely to move directly to a new job and also more likely to move from having no job to the next job. Women have a higher hazard for moving to a no-job state. Given that the population in GLHS is born before 1951, this higher hazard can be explained by women dropping out of the labour market due to marriage or pregnancy. The estimated effects for years of education have high-standard errors. More years of education seem to be negatively associated with moving out of a no-job state when controlling for the effects of time and gender. We do not have a good explanation for this effect. Note that when someone is in the no-job state, this does not imply that this individual is looking for a new job. For example, when people are retired, they are also classified in this state.

Goodness of fit for Model $B$ can be checked visually by comparing model-based survival with Kaplan–Meier survivor curves; see Figure 4. The solid lines are the model-based survival predicted using the baseline data only. The survival is defined for a progressive series of sets of states. The first graph is for survival prior to state 2; that is, survival in state 1. The second graph is for survival prior to state 3; that is, survival in states 1, 7 or 2. The graphs for the higher states are defined analogously. The model-based prediction is derived as the mean of the survivor curves for the 201 individuals in the data. The latent-class structure is taken into account by predicting survival for each individual twice: once assuming the individual is in class 1, once assuming the individual is in class 2. The individual survivor curve is then defined as a weighted sum of the class-specific survival where the weights are the estimated probabilities $π_{k}$ for $k \in {1,2}$ . The dashed lines are the Kaplan–Meier survivor curves (with 95% confidence intervals), where the plus-symbols depict right-censored times.

Figure 3

Transition-specific hazards for women born in 1939–1941 and with less than 12 years of education. Black curves for class 1, and grey curves for class 2 (with 95% confidence bands)

For survival prior to state 2 and for survival prior to state 3, the model-based curves are close to the Kaplan–Meier curves. For the higher states, Figure 4 shows some lack of agreement at the later times part of which is due to the right-censoring. Using the non-parametric Kaplan–Meier estimation as a data summary, the overall similarity with the model-based prediction validates the main features of the fitted model.

Hazards such as the ones in Figure 3, inform about employment history in the sense that they show the differences between the classes and the rates of moving to a next job. Transition probabilities can be used for additional inference. For example, to compare the two classes, we can look at the probabilities for changing states within a 5-year period. Say we do this for a man in birth cohort 1949–1951, with more than 12 years of education, who has been in the labour market for 5 years and is currently in his second job. The probability to move to (or stay in) state $s \in {2,3, 4,5, 6}$ is given by

p_{2 s} (t_{1} = 60, t_{2} = 120 | class = k, (x_{1}, x_{2}) = (0, 1), x_{3} = 1, x_{4} = 1),

Figure 4

Comparison of predicted model-based survival and Kaplan�Meier survivor curves. Solid lines for mean of predicted survival given baseline GLHS data for the 201 individuals. Dashed lines for the Kaplan�Meier curves (with 95% confidence intervals)

which we will denote $p_{2 s}$ , for short. If this man is in class $k = 1$ , then we get

p_{22} = 0.231, p_{23} = 0.343, p_{24} = 0.226, p_{25} = 0.097, and p_{26} = 0.025.

If this man is in class $k = 2$ , then we get

p_{22} = 0.572, p_{23} = 0.283, p_{24} = 0.052, p_{25} = 0.005, and p_{26} = 0.001.

For the computation of these probabilities, we used a piecewise-constant approximation (defined by a 1-month grid) to take the non-constant hazards into account. These transition probabilities illustrate the differences between the classes on the scale of probabilities. For example, we see that for staying in the second job the difference between the classes was already quite large (0.231 vs. 0.572 for classes 1 and 2, respectively). When in class 1, the man has a consistently higher chance to be in another job after 5 years then when he is in class 2, illustrating that class 1 is the class of the movers.

If we want to allocate the above specified man (with identifier i, say) to class $k \in {1,2}$ , we can compute the probability

P (i \in class k | t_{i}, δ_{i}, y_{i}, x_{i}; θ = \hat{θ}) = \frac{g (y_{i n_{i}} ..., y_{i 2} | y_{i 1}, i \in c l a s s k; θ = \hat{θ}) P (i \in c l a s s k | θ = \hat{θ})}{\sum_{ℓ = 1}^{K = 2} g (y_{i n_{i}} ..., y_{i 2} | y_{i 1}, i \in c l a s s ℓ; θ = \hat{θ}) P (i \in c l a s s ℓ | θ = \hat{θ)}},

where $g (y_{i n_{i}} . . ., y_{i 2} | y_{i 1}, i \in c l a s s k; θ = \hat{θ})$ is short-hand notation for the likelihood contribution for individual i conditional on being in class k; see Section 5. For example, if i is a job hopper in the sense that he changes job every year directly with $t_{i} = (0,12,24,36)$ and $y_{i} = (1,2, 3,4)$ , then the probability to be in class 1 is estimated at 0.90 and we would allocate i to class1. However, if i changes job every 5 years, that is, $t_{i} = (0,60,120,180)$ with the same $y_{i}$ , then this probability is estimated at 0.36 and we would allocate i to class 2.

To summarize the GLHS data analysis, Model $A$ shows that inference can be improved by fitting generalized gamma distributions in a model for employment history with time since entry into the labour market as the time scale. Effects of cohort, gender and education are important and have to be taken into account in the statistical inference. In addition, we illustrated that the latent-class Model $B$ is an interesting extension of Model $A$ and allows to explore an assumed latent-class structure. Fitting a two-class structure shows a clear distinction between stayers and movers with respect to employment history. One possible way to extend the analysis would be to investigate covariate effects for class membership. For example, the probability to be a mover ( $π_{1}$ ) could be modelled using logistic regression with gender as a covariate.

7 Conclusion

This article illustrates the use of the generalized gamma distribution for a progressive multistate process in demography. In addition, it is shown that an extension to a latent-class structure is possible and can lead to improved statistical inference.

Because the multistate process in the application is progressive, the likelihood function for the time-dependent model can be constructed using transition-specific hazard functions and state-specific survivor functions. For multistate processes where back-and-forth transitions between the states are possible, the distributional options are the same, but the estimation becomes more complex. For such a case, we would propose to use a piecewise-constant approximation to the continuous-time parametric shape. The approximation would consist of a series of exponential models with changing hazard specifications; see, for example, (Blossfeld and Rohwer, 2002) and Van den Hout (2017).

This article combines parametric hazard models with a discrete distribution for unobserved heterogeneity. Other options are possible. Putter et al. (2007) discuss a semi-parametric alternative for the fixed-effects models, and Hougaard (2000) discusses parametric choices for the distribution for unobserved heterogeneity. Advantages of parametric hazard models are efficiency and the option to predict the process outside the range defined by the data. With respect to the discrete distribution for unobserved heterogeneity, as long as this distribution is defined by a few classes, this option is computationally advantageous and, of course, it circumvents the need to specify a parametric shape for a latent characteristic.

Acknowledgements

The creation of Figure 2 is based on R code on www.everydayanalytics.cawww.everydayanalytics.ca, which is run by Myles Harrison (assessed on 24 September 2017). The authors would like to thank two anonymous reviewers for their detailed feedback.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. Funding The authors received no financial support for the research, authorship, and/or publication of this article.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

References

Aitkin

(1999) A general maximum likelihood analysis of variance components in generalized linear models. Biometrics , 55, 117–28.

Bartolucci

Farcomeni

Pennoni

(2012) Latent Markov Models for Longitudinal Data . London: Chapman & Hall/CRC.

Bijwaard

(2014) Multistate event history analysis with frailty. Demographic Reserach , 30, 1591–620.

Blossfeld

Rohwer

(2002) Techniques of Event History Modeling. New Approaches to Causal Analysis, second edition . Mahwah: Lawrence Erlbaum Associates.

Bongaarts

(2005) Long-range trends in adult mortality: Models and projection methods. Demography , 42, 23–49.

Cole

Bonetti

Zaslavsky

Gelber

(2005) A multistate Markov chain model for longitudinal, categorical quality-of-life data subject to non-ignorable missingness. Statistics in Medicine , 24, 2317–34.

Commenges

(1999) Multi-state models in epidemiology. Lifetime Data Analysis , 5, 315–17.

Cox

Chu

Schneider

Muñoz

(2007) Parametric survival analysis and taxonomy of hazard functions for the generalized gamma distribution. Statistics in Medicine , 26, 4352–74.

Dias

Willekens

(2005) Model-based clustering of sequential data with an application to contraceptive use dynamics. Mathematical Population Studies , 12, 135–57.

10.

Frydman

(1984) Maximum likelihood estimation in the mover-stayer model. Journal of the American Statistical Association , 79, 632–38.

11.

Hougaard

(2000) Analysis of Multivariate Survival Data . New York: Springer.

12.

Jackson

(2011) Multi-state models for panel data: The msm package for R. Journal of Statistical Software , 38, 1–29.

13.

Jackson

(2016) flexsurv: A platform for parametric survival modeling in r. Journal of Statistical Software , 70, 1–33.

14.

Kalbfleisch

Lawless

(1985) The analysis of panel data under a Markov assumption. Journal of the American Statistical Association , 80, 863–71.

15.

Kay

(1986) A Markov model for analysing cancer markers and disease states in survival studies. Biometrics , 42, 855–65.

16.

Mandel

(2013) Simulation-based confidence intervals for functions with complicated derivatives. The American Statistician , 67, 76–81.

17.

Mayer

(2015) The german life history study: An introduction. European Sociological Review , 31, 137–43.

18.

Mueller

Nusbaum

Rose

(1995) The gompertz equation as a predictive tool in demography. Experimental Gerontology , 30, 553–69.

19.

Olshansky

Carnes

(1997) Ever since gompertz. Demography , 34, 1–15.

20.

Prentice

(1974) A log gamma model and its maximum likelihood estimation. Biometrika , 61, 539–44.

21.

Putter

Van Houwelingen

(2015) Frailties in multi-state models: Are they identifiable? Do we need them? Statistical Methods in Medical Research , 24, 675– 92.

22.

Putter

Fiocco

Geskus

(2007) Tutorial in biostatistics: Competing risks and multi-state models. Statistics in Medicine , 26, 2389–430.

23.

R Core Team (2013) R: A Language Environment for Statistical Computing . Vienna Austria: R Foundation for Statistical Computing. URL http://www.R-project.org.

24.

Satten

(1999) Estimating the extent of tracking in interval-censored chain-of-events data. Biometrics , 55, 1228–31.

25.

Stacy

(1962) A generalization of the gamma distribution. The Annals of Mathematical Statistics , 33, 1187–92.

26.

Van den Hout

(2017) Multi-State Survival Models for Interval-Censored Data . Boca Raton, FL: Chapman & Hall/CRC.

27.

Vermunt

Langeheine

Böckenholt

(1999) Discrete-time discrete-state latent Markov models with time-constant and time-varying covariates. Journal of Educational and Behavioral Statistics , 24, 179–207.

28.

Willekens

(2014) Multistate Analysis of Life Histories with R . New York: Springer.

Flexible parametric multistate modelling of employment history

Abstract

Abstract

Keywords

1 Introduction

2. German Life History Survey

Figure 1

Eleven-state process for employment history in GLHS

For the 201 individuals in GLHS, a two-dimensional representation of the start time of the first job and the age at the start of the first job

State table for GLHS: number of times for each successive pair of states. The states are defined by employment history; see Figure 1

3.1 Exponential and Gompertz distribution

3.2 Generalized gamma distribution

5 Estimation

6 Data analysis

Parameter estimates for the mover/stayer model for the GLHS data on employment history (Model B ). Time scale is months since entry into the labour market. Estimated standard errors in parentheses. Estimated intercepts are not included

Transition-specific hazards for women born in 1939–1941 and with less than 12 years of education. Black curves for class 1, and grey curves for class 2 (with 95% confidence bands)

Comparison of predicted model-based survival and Kaplan�Meier survivor curves. Solid lines for mean of predicted survival given baseline GLHS data for the 201 individuals. Dashed lines for the Kaplan�Meier curves (with 95% confidence intervals)

Acknowledgements

Declaration of conflicting interests

Funding

References

Parameter estimates for the mover/stayer model for the GLHS data on employment history (Model $B$ ). Time scale is months since entry into the labour market. Estimated standard errors in parentheses. Estimated intercepts are not included