Healthy longevity from incidence-based models: More kinds of health than stars in the sky

Background. Healthy longevity (HL) is an important measure of the prospects for quality of life in ageing societies. Incidence-based (cf. prevalence-based) models describe transitions among age classes and health stages. Despite the probabilistic nature of those transitions, analyses of healthy longevity have focused persistently on means ("health expectancy"), neglecting variances and higher moments. Objectives. Our goal is a comprehensive methodology to analyse HL in terms of any combination of health stages and age classes, or of transitions among health stages, or of values (e.g., quality of life) associated with health stages or transitions. Methods. We construct multistate Markov chains for individuals classified by age and health stage and use Markov chains with rewards to compute all moments of HL. Results. We present a new and straightforward algorithm to create the multistate reward matrices for occupancy, transitions, or values associated with occupancy or transitions. As an example, we analyse a published model for colorectal cancer. The possible definitions of HL in this simple model outnumber the stars in the visible universe. Our method can analyse any of them; we show four examples: longevity without abnormal cells, cancer-free longevity, and longevity with cancer before or after a critical age. Contribution. Our methods make it possible to analyse any incidence-based model, with any number of health stages, any pattern of transitions, and any kind of values assigned to stages. It is easily computable, requires no simulations, provides all the moments of healthy longevity, and solves the inhomogeneity problem.

Appendix C Some models for the value of fractional occupancy 31 Appendix D Prevalence is a special case of incidence 32 Appendix E One is a special case of many 33 Appendix F Reproduction is special case of value 35 Bibliography 36 1 Introduction 1.1 Longevity and healthy longevity 23 Healthy longevity is an inclusive term for a class of health metrics that combine informa- The matrix U contains probabilities of transition among the transient states; the (i, j) entry The transition matrix P captures the dynamics of the probability that an individual will 170 be in a specified state at t + 1, conditional on its state at t. Let π(t) be a probability vector 171 (non-negative, and entries sum to 1) where π i (t) is the probability that the individual is in 172 state i at time t. Then 173 π(t + 1) = Pπ(t).
(2) 174 In an absorbing Markov chain, eventual absorption (i.e., eventual death) is certain. Before intensity matrix Q has non-positive entries on its main diagonal, and non-negative entries 192 elsewhere, and its column sums are zero.

193
A model with discrete health states and continuous age variation is an inhomogeneous 194 continuous-time Markov chain, for which no solution is available. Thus it is essential to 195 discretise the model for analysis. If Q can be treated as constant over a time interval 196 (t, t + ∆t), then the discrete transition probability matrix over that interval is 198 where e Q is the matrix exponential function (not to be confused with the matrix containing 199 the exponentials of the individual entries). P is called the discrete skeleton of the continuous 200 model by Iosifescu (1980). This formulation requires only the that Q is constant between t 201 and t + ∆t; choosing a sufficiently small time step will make this as true as desired.  4 An alternative commonly used approximation requires, in addition to the assumption of constancy, the assumption that the elements of Q are sufficiently small that a first-order approximation is suitable (Rogers and Ledent, 1976). For transitions over an interval of length ∆t, the approximation is See Schoen (1988) for a comparison of approaches to deriving discrete transition matrices from continuous rate matrices.

Constructing the multistate model 223
The components of the model are a set of transition matrices among health stages with 224 transition probabilities specific to each age class, and an age-advancement matrix that moves 225 surviving individuals to the next older age class. The multistate model is constructed from 226 these components using the vec-permutation matrix formalism, described in detail by Caswell 228 It helps to begin by carefully defining some aspects of the model. 229 τ number of transient (living) health stages α number of absorbing (dead or otherwise removed) stages ω number of age classes s total number of states; s = τ ω + α U x stage transition matrix for age class x, for x = 1, . . . , ω D j age advancement matrix for stage j, for j = 1, . . . , τ M x mortality matrix for age class x, for x = 1, . . . , ω K vec-permutation matrix; parameters τ, ω

230
The matrix U x describes transitions among the health stages for individuals in age class 231 i. The structure of U x is specific to the health condition under consideration; an example 232 will be given in Section 3.2.

233
The age advancement matrix D j moves surviving individuals to the next oldest age class.

234
For example, if ω = 3, then The (ω, ω) entry in the lower right corner is optional, as denoted by the square brackets.

237
If it is 0, then all individuals die after the last age. If it is 1, then the last age class is an 238 open ended interval, with individuals continuing to experience the transitions and mortality 239 defined by U ω . 5

240
The population (or cohort) structure at any time is given by the vector

9
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. matrix, define block-diagonal matrices, of dimension (τ ω × τ ω), Then the block-structured matrix that forms the transient part of the age×stage-classified 248 transition matrix is Here, the matrix K is a special permutation matrix known as the vec-permutation matrix 251 (Henderson and Searle, 1981); see AppendixB for calculation.

252
The mortality matrix M is of dimension α × τ . The (i, j) entry of the matrix M x is the 253 probability of transition from the jth transient state to the ith absorbing state, conditional 254 on age x. The matrixM combining these matrices is The transition matrix for the age-stage classified Markov chain is then   Healthy longevity is an occupancy time and our goal is to compute the statistics of the time 292 spent in some specified set of health conditions. We confront a mathematical problem here. 293 We are, in general, interested in occupancy not of a single state of the multistate Markov 294 chain, but a set of states, defined by certain combinations of ages and health stages. Classical

295
Markov chain theory provides all the moments of the occupancy time for any single state 296 but not for sets of states (e.g., mildly or severely disabled at ages between 60 and 70) . Thus 297 classical theory could provide the statistics of longevity for, e.g., early preclinical CRC at 298 age 60, but not for more interesting combinations such as cancer-free longevity between 60 299 and 70. 300 We solve this problem 6 by using Markov chains with rewards (MCWRs); these permit a 301 stochastic analysis of healthy longevity for arbitrary combinations of ages and health stages. 302 We will see below what this implies in practice for the CRC model, but first we turn to the 303 basic methodology of MCWRs.

11
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

318
The reward r ij is a random variable and is characterised by its moments. Reward matrices 319 can be estimated in many ways, depending on the nature of the reward (van Daalen and 320 Caswell, 2017). For the incidence-based health model, we require rewards for transitions 321 between all age-stage combinations.

322
Since the elements of the reward matrices correspond to the entries of the transition 323 matrix, we write the reward matrices with a block structure corresponding to that of the 324 transition matrix P. We defineR m as the reward matrix containing the mth moments of  Lifetime rewards. Rewards are accumulated from any starting state until absorption.

331
No rewards are accumulated after death (i.e., there is no healthy longevity after death; this 332 seems like a safe assumption). The moments and the statistics (mean, variance, coefficient Theorem 1). The first three moments are, in our notation, whereÑ is the fundamental matrix, The (i, j) entry ofÑ is the mean time spent in state i, before eventual absorption, conditional  The notation here differs from that in van Daalen and Caswell (2017). In that paper, the tilde onρ denoted a vector from which the absorbing states had been removed. Here, however, the tilde denotes vectors and matrices with the block structure of stages within age classes, as given by (6).

12
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 18, 2021. ; https://doi.org/10.1101/2021.04.16.21255628 doi: medRxiv preprint The first three moments suffice to calculate the mean, variance, and skewness of healthy 350 longevity:

356
The vectorρ m contains the mth moments of healthy longevity for all combinations of 357 initial age and health stage. To obtain the moments of healthy longevity as a function of 358 age for a specified initial stage i, we calculate where e i is the ith unit vector of length τ . Similarly, the moments as a function of stage for 361 a specified initial age x are given by where e x is the xth unit vector of length ω.

364
Health and healthy longevity are protean concepts with many definitions. Each has its 365 own reward structure, and in the next section we show how to calculate the componentsB m 366 andC m of the reward matricesR m in equation (12).

367
5 Health reward structures and defining healthy longevity 368 We consider four ways to define healthy longevity.     . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021.  to, you can. Note that total longevity, independent of health status, is given by inserting a 415 1 in every entry of H.

416
8 It is understood that interest might focus on longevity subject to some illness or condition that is not pleasant or desirable.
14 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021. ; The reward indicator vector. The next step is to create an occupancy indicator vector by applying the vec operator to H (stacking the columns on top of each other); This 0-1 vector inherits the block structure of the vectorñ in (6) and the age-stage matrices 420Ũ , etc. Also define a vector ¬h that contains the logical complement of the entries of h; i.e., Partial rewards for partial occupancy. If an individual spends the entire interval within   is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021. ; https://doi.org/10.1101/2021.04.16.21255628 doi: medRxiv preprint first year of life, when infant mortality shifts the fractional occupancy to something smaller, 433 perhaps 0.1 (Chiang, 1968;Preston et al., 2000). 434 We apply the same logic to partial occupancy here. We define a as the fractional occu-435 pancy experienced by an individual that enters or leaves the set S during the interval, and 436 set a = 1 2 .

437
The reward matrices for longevity. We need the moment matricesB m andC m that 438 appear in (12). The calculation ofB 1 must account for three kinds of occupancy: Recalling that h indicates the states in S and ¬h the states in ¬S, the matrix is The first term gives a unit reward to an individual starting and ending the interval in any 445 of the states within S. The second term gives a reward of 1 2 to an individual that starts the 446 interval within S but finishes the interval outside S. The third term gives a reward of 1 2 to 447 an individual that starts the interval outside of S but finishes the interval within S. Because 448 h + ¬h = 1, the result (29) reduces to the (remarkably simple) Transitions into absorbing states from a state in S earn a fractional occupancy reward; 451 the componentC 1 is given by The matrix of first moments of rewards,R 1 , is assembled fromB 1 andC 1 as in (12).

454
Because occupancy is a fixed reward (sensu Caswell 2011; van Daalen and Caswell 2017) the 455 second and third moment matrices are and, in general, This provides all the pieces needed to compute the moments of healthy longevity using  1 that indicate the transitions contained in where i and j denote transient states. Transitions into absorbing states are counted by where i is an absorbing state and j is a transient state.

485
In the case of the CRC model, counting entries into clinical CRC would yield, for any The age-stage matrixB 1 of first moments is obtained using the same vec-permutation 489 construction that producedŨ, and D is as defined in (8). The matrixC 1 is The reward matrixR 1 is assembled as in (12). Because transitions, like occupancy, are fixed 496 rewards, the higher moments are given in equation (34).
497 9 In this model, it is impossible to leave the clinical CRC stages, other than by death. Counting the transitions might be of limited interest, but it is easy to imagine cases where transitions would be more flexible.

17
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  where V m (i, j) is the mth moment of the value of a unit of time in stage i at age j.
Define moment vectors in the same way as the occupancy indicator vector h., Let us assume that partial occupancy for 1/2 a time unit has a value 1/2 of the random Transitions into absorbing states accumulate value only for the fraction of the interval in the 520 state of origin, so that The reward matricesR m are assembled fromB m andC m as in equation (12). Rewards 18 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021. ; https://doi.org/10.1101/2021.04.16.21255628 doi: medRxiv preprint and benefits will generally vary with age, and we would like to assign values to health stage 537 transitions as a function of age.

538
Transitions are, as in Section 5.2, defined by an origin stage, a destination stage, and an 539 age class. The set of transitions of interest at age x is T x . Let v(i, j, x) be the value attached 540 to a transition j → i at age x. The matrices B (x) and C (x) defined in equations (38) and if i and j are transient states, and when i is an absorbing state and j is a transient state.

546
Given these age-specific matrices, the block-structured matricesB m andC m are calcu-547 lated using the vec-permutation matrix as in (41) and (43). The reward matrices are then 548 given by equation (12).

568
The first three moments of remaining longevity are given by the entries of the vectors whereÑ is the fundamental matrix given in (16) andη m is the (block-structured) vector of 573 the mth moments of longevity (Iosifescu, 1980;Caswell, 2013).

574
Suppose, as is true in the example in Section 6, that we distinguish two causes of death, 575 and we are interested in years of life lost to cause 1. In the transition matrix (11), the 576 probabilities of transition to death from cause 1 appear in row 1 of the mortality matrixM.

19
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Following the construction in equation (12) gives the reward matrices The age×stage-structured vectorρ m computed fromŨ,P, and the reward matrices contains 583 the mth moments of life lost, due to death from cause 1, of every age-stage combination. a much longer time in these CRC conditions than will typical individuals, and a focus on 618 mean longevity seriously underestimates these extreme cases.

619
An intuitive feel for this uncertainty is provided by the following calculation. Imagine 620 picking, at random, two individuals identical in age class and health stage. How similar will 621 their eventual healthy longevities be? What is the chance that those longevities will differ 622 by a factor of 2? or a factor of 10? The higher the probability of a k-fold difference, the

22
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021.

25
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  The concept of healthy longevity takes it as given that an individual will live for some length 710 of time and will experience various health conditions during that time. The distinction 711 between incidence-based and prevalence-based analyses reflects the information incorporated.

712
Prevalence-based analyses (the Sullivan method and its relatives) incorporate age, but are 713 necessarily blind to the actual transitions among health stages, relying instead on age-specific 714 prevalences. Incidence-based analyses explicitly include transitions among health stages. The 715 great value of prevalence models is that they can be constructed from cross-sectional data, 716 making them an essential tool for exploring healthy longevity when longitudinal data are not 717 available. The analyses are almost always restricted to expectations (health expectancy). 718 Caswell and Zarulli (2018) presented a replacement of the Sullivan method that provides 719 variances and higher moments using MCWR methods, treating prevalence as an age-specific 720 reward.

721
It is possible, and an interesting exercise, to write prevalence-based models as a degenerate 722 special case of incidence-based analyses. Doing so might bring some extra insight even with 723 limited data; see Appendix D.

724
Health has many dimensions (physical, mental, emotional, spiritual). It may be defined in 725 terms of specific diseases (e.g., colorectal cancer) or in terms of more general conditions (e.g.,  stage (τ = 1) or a single age class (ω = 1); see Appendix E 733 We have emphasized the importance of variance (and higher moments) of healthy longevity.

734
In previous studies, analysis of the variance in longevity has been frustrated by the need to 735 analyse sets of states rather than single states. Classical Markov chain theory provides the 736 moments of occupancy for single states, but not for sets of states. The expected value of a 737 26 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021. ; sum is the sum of the expected values, but the variance of the sum is not the sum of the variances (except in the special case where the items are independent). The MCWR method we introduce here solves this difficulty, not only for occupancy but for transitions and for 740 values associated with them. 11

741
The computational speed of the MCWR calculations will make it easier to use resam-

27
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. is free to proceed in whatever way is most appropriate to the data.

799
A framework for reporting results. There is currently no standard format for reporting 800 multistate models. Every paper seems to do so differently. In our experience most of The analysis presented here invites extensions. We discuss a few here.  Variance components. The variance in healthy longevity in a mixed cohort has two com-823 ponents. One is the variance within each health stage, due to individual stochasticity.

824
The other is the variance between health stages, due to heterogeneity in the rates ex-825 perienced by those individuals. The within-and between-group variance components 826 28 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021. ; quantify the extent to which the variability in outcome is due to within-cohort. See (Caswell, 2014;Hartemink et al., 2017)   If life is a Markov chain, then health is a reward. The combination of this paper and that There are times when a study leads to interesting or entertaining insights that are not a part 864 of the main story but are too good to lose. In this appendix we present some of these.  To our children the matrix methods will seem inevitable because they enable 881 demographers to fall back on a hundred years of development in mathematics -882 due to Sylvester, Frobenius, Kolmogorov and others -rather than to invent the 883 mathematics independently as they go along. 884 Keyfitz may have been overly optimistic, but his points are nonetheless valid. The vec-permutation matrix is, like any permutation matrix, a square matrix with a single 887 1 in each row and column and zeros elsewhere. It rearranges, but does not change the value 888 of, the entries in a vector on which it operates. Let X be an m × n matrix, The vec-permutation matrix K m,n relates the vec of X and the vec of X T :

892
The subscripts on K can be omitted if they are clear from the context. The matrix is given where E ij is a m × n matrix with a 1 in the (i, j) entry and zeros elsewhere (Henderson and 896 Searle, 1981).

897
30 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 18, 2021. ; https://doi.org/10.1101/2021.04.16.21255628 doi: medRxiv preprint C Some models for the value of fractional occupancy 898 In Section 5.3 we presented a model for the reward matrices that specify the value of the 899 partial occupancy of a set of age-stage combinations. We assumed that there was a value for 900 occupancy of an entire interval, with the mth moments given by the vector v m , for m = 1, . . ., 901 and that occupancy for a fixed fraction a of an interval yields a fixed fraction a of the value 902 associated with a full interval (we used a = 1/2).

903
Here, we will derive that result and two other possibly interesting models for assigning . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) 948 and so on. There is no closed form expression for these moments.

949
In terms of the reward matricesR m , this leads to The calculation of health expectancy from prevalence data by the Sullivan method relies on 957 two strong assumptions. The first is that the chance of being healthy at age x, which is given 958 by the prevalence of being healthy at age x, is independent of the health states at x − 1 or 959 x + 1. The second is that the probability of survival from age x to age x + 1 is independent 960 of the health state at x. These assumptions are strong but unavoidable in the absence of 961 longitudinal individual data. It's not the fault of the method.

962
However, independence is a special case of dependence, so it is interesting to deploy our 963 theory in a prevalence-based analysis. Consider the two-state model shown in Figure 10.

964
Let v 1 (x) and v 2 (x) be the prevalences of the healthy and not-healthy states at age x, where 965 v 1 + v 2 = 1, and let p(x) be the survival probability from age x to x + 1. Then the transition 966 matrix U x is 967 U x = p(x)v 1 (x + 1) p(x)v 1 (x + 1) p(x)v 2 (x + 1) p(x)v 2 (x + 1) (80) 968 The equality of the columns of U x means that the health state at x + 1 is independent of the 969 health state at x. The survivors from age x (with probability p(x)) are allocated to healthy 970 and not healthy states according to the prevalences at x + 1. The transition matrix for the 971 32 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 18, 2021.  It is apparent from this formulation how restrictive the hypotheses of prevalence-based 977 analysis are. However, it also suggests some modifications to relax the assumptions. It is 978 simple to extend the number of health stages (e.g., healthy, slightly ill, severely ill, etc.).

979
It is also possible to incorporate differential survival as a function of health state, if such p 1 (x)v 1 (x + 1) p 2 (x)v 1 (x + 1) p 3 (x)v 1 (x + 1) p 1 (x)v 2 (x + 1) p 2 (x)v 2 (x + 1) p 3 (x)v 2 (x + 1) p 1 (x)v 3 (x + 1) p 2 (x)v 3 (x + 1) p 3 (x)v 3 (x + 1) Now the health state at x + 1 is independent of the health state at x, but the survival 985 probability is not. The multistate transition matrixP is constructed as usual, and the array 986 H can now be used to compute healthy longevity in any desired combination of age classes 987 and health stages. 988 We invite researchers confronted with prevalence data to use this analysis as a comple-989 ment to the MCWR analysis of Caswell and Zarulli (2018) and as an alternative to the much 990 more limited Sullivan method.

991
E One is a special case of many 992 We note that one is simply a special case of many; so we can apply the analysis to the very 993 special case of only a single stage or only a single age class. This is as simple as setting 994 τ = 1, or ω = 1, or both. The matrixŨ is, as expected, an age-classified Leslie survival matrix.

1008
The health occupancy matrix H is 1 × ω. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 18, 2021. ; https://doi.org/10.1101/2021.04.16.21255628 doi: medRxiv preprint