Recognizing duration effects in multistate population models

The risk of many demographic events varies by both current state and duration in that state. However, the use of such semi-Markov models has been substantially constrained by data limitations. Here, a new specification of the semi-Markov transition probability matrix in terms of the underlying rates is provided, and a general procedure is developed to estimate semi-Markov probabilities and rates from adjacent population data. Multistate models recognizing marriage and divorce by duration in state are constructed for United States Females, 1995. The results show that recognizing duration in the married and divorced states adds significantly to the model’s analytical value. Extending the constant-α method to semi-Markov models, 2000–2005 U.S. population data and 1995 cross-product ratios are employed to estimate 2000–2005 duration-dependent transfer probabilities and rates. The present analyses provide new relationships between probabilities and rates in semi-Markov models. Extending the constant cross-product ratio estimation approach opens new sources of data and expands the range of data susceptible to state-duration analyses.


Introduction
Multistate models typically follow persons as they move from state to state over age and/or time, using the Markov assumption that the risk of movement depends only on a person's current status. However, demographers have long known that, in many situations, the duration or length of time a person has been in their current state, can substantially affect the risk of an interstate movement.
Models that recognize both current state and duration in that state are known as semi-Markov models, and a substantial body of statistical, actuarial, and demographic literature has explored them. Life insurance actuaries have long used "Select and Ultimate" life tables, which are based on mortality by age, sex, and years since the policy was purchased. In the USA, an ongoing Society of Actuaries mortality investigation used a 15-year select period (Jordan Jr., 1975, p. 24-28). Actuarial experience showed that individuals at short policy durations had higher death rates than similar persons at long durations. Schoen (1977) showed that the probability of divorce varied with both age and duration of marriage, with duration having a greater effect than age.
The statistical properties of semi-Markov models were explored in depth beginning in the 1950s (e.g., Smith 1955), and that work has since been extended (e.g., Cinlar, 1969;Feller 1968;Ginsberg, 1971;Grabski, 2016;and Hoem, 1972). With regard to multistate models, the pioneering work of Wolf (1988) was the first to examine life table construction with duration dependence. Other noteworthy contributions/applications were made by Cai et al., 2006;Cook and Lawless (2018); Hennessey (1980); Keilman and Gill (1986); Lynch and Brown (2010); and Rajulton (1985). The Appendix section provides additional references and citations to implementing software.
Here, there are two objectives. First, to relate multistate semi-Markov rates to semi-Markov probabilities and explore the implications of those relationships. Second, to extend the Constant-α approach to estimate semi-Markov probabilities and rates from data on adjacent populations in the context of a never married/married/divorced multistate model that recognizes duration of marriage and duration since divorce.

Specifying the state-duration transition probability matrix
To develop the semi-Markov model, let there be S states, where state S i has j durations. Each of the N state-duration categories is termed a cell. For simplicity, assume no mortality or other source of attrition, and age/time intervals of n years. To be concrete, and with little loss of generality, we emphasize the 2-state model where each state has 3 specified durations: 0, 1, and 2, the first two of n years, and the last duration category open-ended.
The logical structure of an age-duration model means that the possible flows between cells are restricted. Persons at each duration in each state can only move (that is move between states, not advance to a higher duration) to duration category zero in one of the other S−1 states. Accordingly, in the 2-state, 3-duration, model of Eq. (1), there are 2 (origin states) × 3 (durations) × 1 (destination state) = 6 possible transition rates.
Markov models are generally based on the underlying forces of transition or the corresponding rates of transfer from one state to another. In the semi-Markov context, however, persons also advance to a higher duration in the same state, a status change that is not directly captured by any occurrence/exposure transfer rate. Semi-Markov probabilities reflect both interstate transfer and intrastate advancement. Specifically, with N = 6, the 6 × 6 transition probability matrix of a 2-state, 3-duration model has the form A ¼ π 10;10 π 11;10 π 12;10 π 20;10 π 21;10 π 22;10 π 10;11 0 0 0 0 0 0 π 11;12 π 12;12 0 0 0 π 10;20 π 11;20 π 12;20 π 20;20 π 21;20 π 22;20 0 0 0 π 20;21 0 0 0 0 0 0 π 21;22 π 22;22 2 6 6 6 6 6 6 4 3 7 7 7 7 7 7 5 where π hi,jk is the probability that a person in state h at duration i at the start of an interval is in state j at duration k at the end of the interval. The rows of A represent destination states, while the columns of A represent origin states. In the model we are considering, where there is no attrition, all columns sum to 1. Of the 36 cells of matrix A, only 18 have non-zero probabilities. There are 3 non-zero probabilities in each column, as each person can, at the end of the interval, either advance to the next higher duration (or stay at the highest duration) in the initial state or be at duration zero in any state. Multiple moves (i.e., changes of state) in an interval are possible, and a person can be at duration zero in the initial state by moving to another state and then returning. All persons who move during an interval are at duration zero at the end of the interval.
In much life table construction, the transition probability matrix is found from a matrix of underlying rates. Because a conventional rate matrix does not capture advancement, going from transfer rates to semi-Markov probabilities requires a cell by cell rather than a matrix-to-matrix approach. Nonetheless, given a set of occurrence/ exposure transfer rates, every probability in A can be expressed in terms of those rates using established multistate calculation procedures.
To do so, we follow the procedures in Schoen (Schoen, 1988, Chap.4;Schoen, 2006, Chap.1). It is convenient to start with the (2,1) cell of A, where π 10,11 is the probability that a person initially in state 1 at duration 0 (i.e., years 0 through 4) will be in state 1 at duration 1 at the end of the interval. That probability of advancement to duration category 1 is simply the probability of never leaving state 1, or where m hi,j0 is the rate of movement from state h, duration i to state j, duration 0, and the implicit survivorship function is assumed to be linear (see also Preston et al., 2001, Chap. 3). This approach serves to specify all of the six π's in rows 2, 3, 5, and 6 of A. For present purposes, we assume linear survivorship as that assumption is generally reasonable and yields explicit algebraic solutions. Alternative solution procedures, such as those discussed in Schoen (1988), Chap. 4), are possible but would make little substantive difference here. Should a large rate be encountered that makes the linear assumption problematic when 5-year intervals are used, interval length can be reduced to 1 year.
Continuing down the first column of A, π 10,20 is the probability that a person initially in state 1 at duration 0 ends the interval in state 2 at duration 0. To find that semi-Markov probability, we need a multistate calculation. Consider a 2-state Markov model that does not recognize duration but that has transfer rates from the initial stateduration equal to those prevailing in the semi-Markov model. The other transfer rates in this 2-state Markov model are considered equal to those at duration 0 in the semi-Markov model. Then, π 10,20 is the same as p 12 , the Markov model probability of starting in state 1 and ending in state 2. Under the linear assumption, in a 2-state model where m ij is the occurrence/exposure rate of transfer from state i to state j, the matrix of rates is of the form and the transition probability matrix is where p ij is the probability that a person initially in state i will be in state j at the end of the interval. As there is no attrition, each column of M sums to zero, with the diagonal element equal to minus the sum of the off-diagonal rates. Each column of PP sums to one.
Under the linear assumption (Schoen, 2006, Section 1.6), Markov transition probability matrix PP can be found from rate matrix M by where I is the N × N identity matrix (which has ones on the main diagonal and zeros elsewhere). Thus π 10,20 follows from p 12 using rates from duration zero to duration zero. The algebraic expression is π 10;20 ¼ 2n m 10;20 = 2 þ n m 10;20 þ m 20;10 The values of π 11,20 and π 12,20 can also be obtained from the p 12 element, but by using rates from state 1 at durations 1 and 2, respectively. All first transfer rates take persons to the other state at duration zero, while subsequent rates take those persons from duration zero to duration zero.
To find π 10,10 , the remaining probability in the first column, we need only use the fact that the sum of each column is one, hence π 10;10 ¼ 1−π 10;11 −π 10;20 ð7Þ Algebraically, the result is π 10;10 ¼ 2n 2 m 10;20 m 20;10 À Á = 2 þ n m 10;20 À Á 2 þ n m 10;20 þ n m 20;10 The same approach yields all of the remaining probabilities in A in terms of the set of underlying rates. By straightforward extension, it takes any set of transition rates and provides the elements of the associated semi-Markov transition probability matrix. Table 1 provides the complete algebraic solution for the probability matrix of a 2-state, 3-duration model in terms of the underlying transition rates.
While Table 1 gives the 18 probabilities in A in terms of the underlying 6 rates, we can also take those 18 equations and solve for the 6 rates and for 12 probabilities in terms of the other 6 probabilities. In the case of probability matrix A of Eq. (1), Table 2 gives explicit solutions for the 6 rates and for 12 probabilities in terms of the other 6 πs. The implications of doing so are substantial and give rise to what may be termed a rate principle: the number of independent probabilities in a semi-Markov model is given by the number of independent rates underlying that model. Matrix A, like its underlying rate matrix, has only 6 independent elements. It follows that an arbitrary probability matrix of the form of A is likely to have no underlying set of semi-Markovian rates because its elements would not be constrained by the relationships in Table 2. In other words, that probability matrix is not "embedded" in a Markovian process (cf. Singer & Spilerman, 1976). The problem of finding rates from such probabilities is addressed in a later section.
Here, let all rates from every given state to another given state be the same at all durations. In the context of our 2-state, 3-duration model, there are only 2 distinct rates, m 12 and m 21 . Using the approach of the preceding section, we can write all of the elements of the semi-Markov transition probability matrix in terms of those two rates. The dominant right eigenvector of that probability matrix provides the long-term (stable population) state-duration composition (Schoen, 2006). That eigenvector, u, can readily be found from A using mathematical software such as Maple or Mathematica.

Estimating transition probabilities from adjacent populations under constant-α
There are a number of situations where population figures by state and duration at both the beginning and end of an age/time interval are known, but there is no information on the transitions during that interval. The constant-α approach, presented in Schoen (2020), can be extended to the semi-Markov case and allow interstate probabilities to be estimated. This section describes how to do so.
The constant-α approach is based on the assumption that the cross-product ratios (α's) of the multistate transition probability matrix are fixed. Cross-product ratios are analogous to odds ratios, can be formed from any rectangular set of 4 non-zero matrix elements, and equal the product of the upper left and lower right elements divided by the product of the lower left and upper right elements. For example, in A, we can define α 1142 ¼ π 10;10 π 11;20 = π 10;20 π 11;10 À Á which is one of the 7 distinct α's in A. The subscripts "1142" represent the upper-left (1,1) and lower right (4,2) elements of the ratio. A distinct cross-product ratio includes at least one cell that is not included in any other cross-product ratio.
If the transition probability matrix is viewed as a contingency table, the constant α's can be interpreted as the fixed interaction effects of a saturated log linear model. Preserving α's can provide maximum likelihood estimates that maximize entropy, as they find the pattern of interstate flows that can arise in the greatest number of ways. In multistate Markov models, Schoen (2020) described how to estimate transition probabilities from a variety of data sources and found that the approach provided good estimates of movements between poverty states in the USA.
Here, we seek to implement the constant-α approach in the semi-Markov context where data are available on adjacent populations. Let x jk represent the start of interval population in state j at duration k, and let y jk represent the end of interval population in state j at duration k. Then where P, which has the form of A, is the transition probability matrix and vectors x and y contain the x jk and y jk population values, respectively.
In the no-mortality semi-Markov case, let us rewrite Eq. (11), using a base transition probability matrix, B, whose elements imply the set of cross-product ratios that are being held constant. Matrix B should be chosen with care and needs to reflect a population with the same state-duration structure and the same interstate movements as the population whose probabilities are to be estimated.
To satisfy the projection relationship, matrix B is pre-multiplied by a diagonal matrix R of row factors and post-multiplied by a diagonal matrix, C, of column factors. The ith diagonal element of R is r i, with r 1 = 1, and the j-th diagonal element of C is c j . Hence, we can write where the desired transition probability matrix, P, is given by where z can be either x or y. By the definition of α, matrix P has the same crossproduct ratios as matrix B. However, the elements of P generally do not satisfy the constraints of Table 2 even when the elements of B do.
With N state-durations, Eq. (12) has (2N − 1) unknowns, the N diagonal elements of C and (N − 1) diagonal elements of R. Those (2N − 1) unknowns can be found from the (N − 1) independent scalar projection equations contained in Eq. (12) and the N equations that require that the N columns of P sum to 1. An iterative solution can be found, but here we proceed by solving the (2N − 1) equations. That approach has the advantage of finding all of the possible solutions. There can be more than one valid (i.e., real and non-negative) demographic solution, while there may be no valid solutions at all. The latter can arise if the cross-product ratios are incompatible with the given populations, the most obvious case being when a large ending population at one duration arises solely from a small initial population at the previous duration.
When the probabilities and rates are known, they can be used in life tables or other demographic models. For example, the life course of a cohort can be traced by a multistate life table, and all of the life table summary measures calculated. We now turn to applying the approaches presented here, first to use rates to calculate a state-duration life table, and second to estimate interstate transfer probabilities and rates using the constant-α method.

Calculating a state-duration model from duration-specific rates
Here, we calculate a semi-Markov model by starting with a Markovian multistate model and extending it through the introduction of duration-specific rates. Marital status models are particularly appropriate for such extensions, as both divorce and remarriage after divorce are known to vary by duration in state.
We begin with the age-state-specific rates used in the construction of the marital status life table for United States Females, 1995 (cf. Schoen & Standish, 2001). To simplify matters, the semi-Markov calculations proceed from age 15 to age 50, ignoring mortality. That yields a 3-state model with states never married (s), married (m), and divorced (v).
We extend the 1995 life table by adding 5-year duration categories 0 and 1, and open-ended duration category 2, to both the married and divorced states. Data on second marriages by duration of first divorce and age at divorce are available for 1995 from Bramlett and Mosher (Bramlett & Mosher, 2001, Table 7) and provide the basis for allocating age-specific remarriage rates (m vm ) to the three duration categories. Ageduration-specific divorce rates (m mv ) for first marriages in California, 1969, are provided in Schoen (Schoen, 1975 Table 2). While somewhat old, they appear to be the most suitable values available. The relative sizes of those published duration-specific rates were then weighted by the initial state composition at each age interval in the extended life table. The weighted differential values, by duration, were multiplicatively adjusted to reproduce the all-durations rate in the 1995 life table. Those adjusted duration-specific rates were the inputs used to calculate the extended multistate life table.
The construction of the extended life table proceeded age by age, beginning with 100,000 persons in the never married state at exact age 15. The state-duration composition of the extended table at the end of each age interval is generated from the initial state-duration composition survived, per Eq. (11), by a 7 × 7 state-duration transition probability matrix. That transition matrix is the 6 × 6 matrix of Eq. (1), with a top row and left-most column added to reflect the never married state. The expressions for the marriage and divorce cells of the matrix are shown in Table 1. There is no re-entry to the never married state, and the probabilities that a never married person ends the interval never married (π ss ), married at duration 0 (π s,m0 ), and divorced at duration 0 (π s,v0 ) are The linear assumption is used throughout. Persons moving between states always begin the next interval at duration 0. The extended life table terminates at exact age 50, after which mortality is more salient and there are fewer marital status transitions. The source 1995 rates and the extended life table functions are given in Table 3 Those results may seem inconsistent at first, but the figures in Table 4, panel B and the first panel of Table 3 offer an explanation. Divorces are rather evenly distributed over the three duration categories, but remarriages are heavily (71%) concentrated at duration 0. Divorce rates decline gradually over age, while remarriage rates drop sharply after age 35. Thus, the 3-duration extended life table has faster and earlier remarriage, which shortens the average duration of a divorce and lengthens the average duration of a marriage. Recognizing duration does make a difference.

Estimating probabilities from adjacent populations using constant-α
The approach here uses the cross-product ratios from the 1995 extended life table of  the previous section to estimate duration-specific probabilities from marital status life  table populations Schoen (2016). Following the procedure described after the presentation of Eqs. (11)-(13), the (2N − 1) = 13 equations were solved for the row and column adjustment factors to the 1995 base probabilities. There were multiple solutions, but only one was demographically appropriate (i.e., with all rates between 0 and 1; though rates can exceed one, such a rate would be unrealistic here). All of the adjustment factors were fairly close to 1, varying only from 0.70 to 1.61. The 2000-2005 estimated matrix of probabilities, P, for ages 30 to 35 was then calculated using Eq. (13). The result is with all columns summing to 1. The largest interstate movement probabilities are from the divorced states to state m 0 . Married persons have probabilities of remaining married of greater than 80%.
In sum, the calculation of an estimated transition probability matrix from a base probability matrix and adjacent populations is straightforward. However, the calculation of the interstate movement rates (and decrements) from the adjacent populations and matrix P probabilities is more complicated and is examined next.

Calculating the non-Markovian marriage and divorce rates and decrements
Estimated transition probability matrix P is non-Markovian because constraints such as those given in Table 2 generally do not hold. Finding appropriate rates consistent with the input populations and estimated probabilities is a non-trivial problem that, to the best of my knowledge, has not been carefully examined in the demographic literature.
In order to find occurrence/exposure rates satisfying Eqs. (12) and (16), more than 7 distinct rates are needed, and there is no unique solution. Here, a 2-step approach is proposed.
Step 1 distinguishes between rates that describe a person's first interstate movement and those that relate to a subsequent movement. Let Mf denote a first move rate, and M a subsequent move rate. To introduce decrements, let df jk be the number of first moves from persons in state-duration j at the start of the interval who move to state k during the interval.
There are 7 first decrement rates, one from each state-duration, and every first move has to be to duration zero in the other state. These rates are related to the probability of a first decrement, and in the linear case can be described by an expression like Eq.
(2) to solve for Mf in terms of π yields where h is the state-duration where persons initially in state-duration j would be at the end of the interval, absent a move. Eq. (17) provides all 7 first transfer rates. Again using established linear life table relationships, the 7 first decrements produced by those rates are of the form where x j is the beginning of interval population in the initial state-duration.
To find the subsequent rates and decrements, it is helpful to set out the 7 stateduration model algebraically by writing 7 equations that describe all of the interstate flows. Those 7 flow equations are The first five flow equations follow from the first decrements as defined above, that is the first movements based on the person's initial state-duration. The move of a person initially in state-duration m 0 who advances to state-duration m 1 and then moves to state-duration v 0 during the interval is included in df m0,v0 , and hence in Mf m0,v0 . Since there is no attrition, summing all of the seven flow equations confirms that the total ending population equals the total initial population. Thus, there are only six independent flow equations.
The last two flow equations are conceptually different and include subsequent moves between state-durations m 0 and v 0 . Those two equations do not include terms for x m0 and x v0 because those persons, absent a move, would be in state-durations m 1 and v 1 , respectively at the end of the interval. All subsequent moves from state-durations m 0 and v 0 must come from entrants during the interval, i.e., the df terms in those flow equations. Under the linear assumption, those entries are, on average, at mid-interval. It follows that (n/2) times the ending (y m0 or y v0 ) population reflects the number of person-years lived in state-duration m 0 or v 0 during the interval. Multiplying those person-years by the M m0v0 or M v0m0 rate of subsequent movement provides the number of subsequent moves between state-durations m 0 and v 0 .
In general, the first and subsequent rates for the same transition differ. Assuming M = Mf produces values that do not satisfy the flow equations. Furthermore, those last two flow equations reveal a further difficulty: they only determine net subsequent decrements, that is the difference [(n/2) у v0 M v0,m0 − (n/2) у m0 M m0,v0 ].
To surmount that difficulty and calculate the subsequent rates and decrements, we go to Step 2. Borrowing from Schoen and Jonsson (2003), we assume that the product of the rates of divorce and remarriage remains constant. The heuristic argument is one of "attractiveness": if (re)marriage becomes more (or less) attractive, one of the rates is likely to rise and the other to fall, so their product can remain unchanged. Thus, we can write Using Eq. (20) with one of the last two flow equations in Eq. (19) allows the calculation of the two subsequent (M) rates and decrements.
The results of the 2-step calculations for the rates and decrements values are shown in Table 5, along with the beginning and ending populations by state-duration. First move divorces occur in roughly equal numbers in the three duration groups, while first move remarriages are concentrated at duration zero. Table 6 summarizes the seven state-duration model. At ages 30 to 35, the cohort of 100,000 women have a total of 13,471 marriages and 8726 divorces. First decrement divorces were 76% of all divorces, while first decrement remarriages were only 59% of all remarriages, a reflection of the high remarriage rates in the years immediately following a divorce.
The 2-step approach presented in this section permits the calculation of rates and decrements from estimated non-Markovian transition probability matrices, such as the one in Eq. (16). While the solution is not unique because there is insufficient information to fully identify the model's non-Markovian aspects, a reasonable, demographically sound solution is presented. These procedures extend the constant-α approach to fully provide semi-Markov probabilities, rates, and decrements from a base probability matrix and adjacent population values.

Summary and conclusion
Semi-Markov multistate models, which recognize both current state and duration in that state, are frequently useful in demographic analyses. The risk of many vital and health events, such as marriage, divorce, and recovery from disability, can vary greatly by duration in state, and that differential risk is often worth examining. Note: Model states are never married (s), married 0-5 years (m 0 ), married 5-10 years (m 1 ), Married 10 or more years (m 2 ), divorced 0-5 years (v 0 ), divorced 5-10 years (v 1 ), and divorced 10 or more years (v 2 ) Source: See text provides a brief, non-technical introduction from an actuarial perspective. Hoem (1972) and Cook and Lawless (2018) are more statistical, but provide good introductory treatments. More advanced treatments can be found in Cai et al. (2006) and Barbu et al. (2017). The computer programs in this paper were written using Maple software, and other mathematical packages, such as Mathematica, can also be used. The computer package R has the most developed semi-Markov software. Willekens and Putter (2014) give an excellent discussion of multistate software in general, with some useful information for semi-Markov modeling. Some specific semi-Markov packages in R are examined in Alvares et al. (2018) and in Krol and Saint-Pierre (2015).