Computation of Multi-State Models using GAUSS, A Matrix Based Programming Language

NOT FOR QUOTATION WITHOUT THE PERMISSION OF THE AUTHORS

Computation of Multi-state Models Using

_GAUSS.

A

Matrix Based Programming Language

Andrew Foster Nut han Ke.JLfitz

December 1986 WP-86-75

Working P a p e r s are interim r e p o r t s on work of t h e International Institute f o r Applied Systems Analysis a n d have r e c e i v e d only limited review. Views or opinions e x p r e s s e d h e r e i n d o not necessarily r e p r e s e n t t h o s e of t h e Institute or of i t s National Member Organizations.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS A-2361 Laxenburg, Austria

Foreword

A g r e a t v a r i e t y of software e x i s t s to help demographers to perform t h e i r cal- culations. What e x c u s e , then, i s t h e r e f o r t h e f u r t h e r set of software t h a t t h i s working p a p e r p r e s e n t s ?

The answer lies in t h e exceptional qualities of GAUSS, a programming language developed f o r micro computers in t h e last y e a r or t w o by Lee Edlefson a n d his associates. I t handles numbers as a r r a y s . So d o o t h e r languages, f o r in- s t a n c e APL. But GAUSS i s simpler, f a s t e r , and more flexible than t h e s e competi- tors. I t uses t h e s t a n d a r d keyboard, a d d r e s s e s a memory of 400k or more, and i s e a s y to l e a r n .

GAUSS i s also v e r y f a s t , taking full advantage of t h e c o p r o c e s s o r of a n XT or a n AT. Speed c a n b e important in multi-state calculations, and in simulations.

GAUSS, f o r instance, creates a n 8 0 x 8 0 matrix of random elements, i n v e r t s i t , and verifies t h e inversion to 16 decimal places, all in less than 1 minute.

Beyond t h e s e advantages i s t h e brevity of i t s programs. The e n t i r e multi- state life table c a n b e written on 25 lines; i t s e s s e n c e a p p e a r s on t h e bottom half of t h e p a g e containing Appendix I. The usual life t a b l e f o r mortality i s a s p e c i a l case to which t h e p r o g r a m i s applicable; i t i s equally applicable to single y e a r s of a g e up to 90. W e have available diskette versions of t h e programs, b u t t h e y are s h o r t enough t h a t they c a n b e keyed in not much more t h a n t h e time r e q u i r e d t o write a letter asking f o r t h e diskette.

On a personal note: a f t e r programming in FORTRAN and in BASIC f o r most of 25 y e a r s , I have now t u r n e d to GAUSS f o r all of my s e r i o u s work. One c a n l e a r n enough t o g e t s t a r t e d in i t within a day; to l e a r n to use i t well and flexibly of c o u r s e t a k e s longer. I am g r a t e f u l to Andrew F o s t e r f o r instructions on GAUSS, and f o r t h e elegance t h a t t h e programs of this working p a p e r exhibit.

Nathan Keyfitz L e a d e r

Population Program

-

ⁱⁱⁱ

-

Abstract

In t h i s p a p e r w e p r e s e n t a s e r i e s of programs in GAUSS. Three areas of ap- plication are shown: (1) construction of a multi-state life t a b l e , (2) projection of a multi-state population, and (3) t h e calculation of v a r i a n c e s and covariances f o r a p r o j e c t e d population. Illustrative r e s u l t s a r e p r e s e n t e d f o r marital s t a t u s d a t a by a g e f o r Canadian women in 1981. Program code and t h e a p p r o p r i a t e d a t a are pro- vided so t h a t t h e r e a d e r may r e p r o d u c e o u r r e s u l t s and develop a working knowledge of multi-state analysis. The programs are readily applicable to d a t a f o r o t h e r countries, and to multi-region as w e l l as to marital s t a t u s analysis.

Contents

Introduction

The Multi-state Life Table Multi-state Projection

Variances in the Multi-state Model References

Appendix I Appendix I1 Appendix I11 Appendix IV Appendix V

Page

Computation of Multi-state Models

Uaing

GAUSS.

A Hatrix Based Programming

Language Andrew Foster*, Nathan K e f l t z * *

Introduction

Despite t h e widespread recognition of t h e r e l e v a n c e of t h e multi-state model in t h e analysis of demographic p r o c e s s e s , t h e application of multi-state methods h a s been limited by t h e s i z e and complexity of computer programs seemingly re- quired f o r such analysis. R e s e a r c h e r s i n t e r e s t e d in gaining a p r a c t i c a l under- standing of t h e s e methods have often been faced with a need to invest a substantial amount of time to develop and modify t h e a p p r o p r i a t e software. In t h i s p a p e r w e p r e s e n t new and simple programs.

The programs are all written in GAUSS version 1.40, a matrix based language f o r t h e IBM or IBM-compatible p e r s o n a l computer. The advantage of a matrix based language i s s e e n in t h e code p r e s e n t e d here-programs are s h o r t and closely re- lated to t h e analytic e x p r e s s i o n s on which they are based. These programs could b e a d a p t e d to o t h e r matrix based languages such as APL or PROC MATRIX in SAS.

The programs are designed to t a k e advantage of t h e PROC command, which al- l o w s one in e f f e c t to design a specialized language by loading into memory previ- ously compiled p r o c e d u r e s . A v a r i e t y of different applications programs may then b e developed which c a l l t h e same g e n e r a l procedures.

Some g e n e r a l i t y and efficiency in t h e p r o c e d u r e s h a s been s a c r i f i c e d in o r d e r to improve readability. W e h a v e avoided t r i c k s t h a t produce d e s i r e d r e s u l t s quick- ly b u t in o b s c u r e ways. In addition, w e have assumed t h a t t h e programs are to b e applied to problems of relatively l o w dimension. GAUSS i s at i t s b e s t when process- ing matrices, b u t t h e maximum matrix s i z e i s at p r e s e n t 6 4 K or about BOO0 ele- ments. This means, f o r example, t h a t a population projection by a Leslie matrix t h a t pre-multiplies a v e c t o r of single-year a g e counts c a n g o only up to a g e 90. An a l t e r n a t i v e method c a n b e c o n s t r u c t e d t h a t does not r e q u i r e t h e full projection

*Andrew F o s t e r , Graduate Croup i n Demography, U n i v e r s i t y o f California, 2234 Piedmont Avenue, B e r k e l e y , CA 94720, USA

**Nathan K e y f l t e , Leader, Population Program, IIASA, A-2361 Laxenburg, A u s t r i a

matrix which is r a t h e r s p a r s e . The simplicity of t h e code, however, would b e s a - c r i f i c e d .

Finally, an attempt h a s been made to permit even f i r s t time u s e r s of GAUSS to r e p r o d u c e some of t h e published results-one need only e n t e r t h e programs (in- cluding t h e d a t a formatting program) in ASCII files with t h e specified names, and t h e n r u n t h e program "msstart.can" from GAUSS. A t t h a t point any of t h e applica- tions p r o c e d u r e s c a n b e run. S e e Appendix V f o r t h e full set of commands.

The multi-state life table

The program which p r o d u c e s a multi-state life t a b l e r e l i e s on t h e development in R o g e r s (1975). so r e a d e r s who wish a more detailed discussion are r e f e r r e d t h e r e ; in t h i s p a p e r w e provide a s h o r t review.

The s t a r t i n g point f o r a multi-state analysis i s t h e construction of a n ap- p r o p r i a t e schedule of transition intensity matrices. If t h e r e are k s t a t e s , t h e n f o r e a c h a g e z w e would h a v e a k X k matrix p ( z ) . The i j ( i

<>

^j ) element of p ( z ) i s minus t h e rate of transition from t h e jth to t h e ith state at e x a c t a g e z . The i i t h element is t h e rate of movement o u t of state i and i s thus t h e mortality rate in state i plus t h e sum of t r a n s i t i o n rates o u t of state i to t h e o t h e r m -1 states at t h a t age. Note t h a t in t h i s formulation d e a t h is not included explicitly as one of t h e states. An a l t e r n a t i v e formulation would include d e a t h as a n absorbing s t a t e , thus increasing t h e dimension of t h e matrices p ( z ) by one.

A s i s t h e case i n t h e simple life t a b l e , t h e multi-state life t a b l e does not work d i r e c t l y with t h e underlying t r a n s i t i o n intensities (supposed continuous with age) but with a n approximation in finite age g r o u p s a n d assuming a constant transition intensity matrix

M,

within a g e g r o u p x . The multi-state life t a b l e p r o c e d u r e , then, t a k e s as input a matrix VCM which is a v e r t i c a l concatenation (v.c.) of t h e s e d i s c r e t e t r a n s i t i o n intensities

4.

A f u r t h e r approximation i s made when w e c o n s t r u c t t h e matrix p,, t h e pro- portion surviving by state from t h e beginning of e a c h a g e g r o u p , z , to t h e begin- ning of t h e n e x t a g e g r o u p , z +h ,

where h i s t h e s i z e of t h e a g e group, t a k e n in t h e examples considered h e r e to be 5 y e a r s . An e x a c t e x p r e s s i o n would r e q u i r e a s p e c t r a l decomposition of t h e matrix -h M,, o r some o t h e r complication. While GAUSS provides a mechanism f o r t h a t

decomposition in t h e form of a g e n e r a l eigenvalue p r o c e d u r e , as long as t h e transi- tion intensities and a g e intervals a r e s m a l l t h e i n c r e a s e in computation time i s dif- ficult to justify.

The calculations used f o r o t h e r life table measures a r e straightforward gen- eralizations of those used f o r constructing simple life tables. Indeed, since t h e

GAUSS

language o p e r a t e s on matrices t h e individual lines look r a t h e r like t h e lines t h a t might a p p e a r in a PASCAL program, f o r example, t h a t i s used f o r t h e single s t a t e case. The matrices L,, survivorship from t h e beginning t o a g e z, a r e con- s t r u c t e d by s t a r t i n g with a n L o (if t h e radix a g e i s 0) which is t h e identity matrix instead of t h e usual s c a l a r 1 and multiplying by successive p, matrices. The person-years in t h e a g e g r o u p hL, are estimated by a l i n e a r approximation within e a c h a g e group. The T, a r e t h e sum of t h e h L y f o r y 2 z. Finally, t h e s t a t e depen- dent life expectancy e , whose i , _{j t h} element i s t h e expected number of y e a r s t o be s p e n t in s t a t e i f o r a n individual c u r r e n t l y in s t a t e j , i s T, X &-I. The matrices VCL, VCLL, VCT, and VCE a r e t h e v e r t i c a l concatinations of t h e LZphL, T,, and e, matrices, respectively.

One version of t h e p r o c e d u r e a p p e a r s in Appendix I. One might readily use a cubic f o r integration to estimate person y e a r s and provide a more r e a l i s t i c closing of t h e life table. Whatever t h e methods f o r integration and closing, given t h e ma- t r i x of d i s c r e t e transition intensities and a radix a g e t h e resulting multi-state life table will not be v e r y d i f f e r e n t from ours.

Appendix I a l s o contains a n applications program f o r t h e durability of m a r - r i a g e in Canada. In o r d e r to facilitate t h e application i t i s helpful t o c r e a t e anoth- e r p r o c e d u r e , MSMASK. MSMASK allows one t o conveniently perform counterfactu- a1 experiments by setting various elements of t h e transition matrix t o z e r o using a f i l t e r o r "mask". Multiplication of the d a t a matrix by MSMASK is element-by- element, not t h e usual row-by-column multiplication. One c a n give interesting sub- stantive i n t e r p r e t a t i o n s t o t h e r e s u l t s obtained when c e r t a i n masks a r e applied.

F o r example, t h e expected time in t h e married s t a t e f o r a married individual in t h e a b s e n c e of r e t u r n s t o marriage from divorce o r widowhood will be t h e expected time s p e n t in the current marriage. In t h e s e programs we consider f o u r types of masks in addition t o 'Base", t h e mask of ones which has no effect: "NoReMa", "No- Div", "NoWid", and "NoDea" to a b s t r a c t from remarriage, divorce, widowhood and death, respectively.

Table 1 provides a n illustration of this example. For two a g e s (20 and 65 with radix 20) we p r e s e n t 1, and e, f o r a base r u n as well as a f t e r abstracting from remarriage. Consider i and j E Is,M,w,DJ (e.g. single, married, widowed, and di- vorced), where j indexes columns o r "originating status" and i indexes rows o r

"end status". If x i s a g e in Table 1 , then then lij(x) may be i n t e r p r e t e d as t h e probability t h a t a n individual originating in state j at age 20 will be in state i _at a g e x ; ei,(x) i s t h e expected number of y e a r s an individual in state j at age x will spend in state i b e f o r e dying. Thus, a woman who is single at a g e 20 may expect 38.7 y e a r s of married life. If w e do not allow remarriage, however, s h e can expect only 32.2 y e a r s of married life. If mortality rates were t h e same in t h e t h r e e ever-married s t a t e s , then t h e a v e r a g e length of time spent in a f i r s t marriage f o r a single woman of a g e 20 is 32.2 y e a r s .

Multi-state projection

One extension of t h e multi-state life table i s the multi-state projection which uses a block-Leslie matrix, where t h e t e r m "block" r e f e r s t o t h e f a c t t h a t each non-zero element in t h e usual Leslie matrix is taken t o be a k ^X k matrix when t h e r e are k states. W e constructed such a block-Leslie matrix f o r Canadian fe- males (Feeney, 1970) and use the matrix in the standard fashion t o provide projec- tions based on r a t e s and counts from 1981. Since d a t a on fertility by a g e and s t a t e were not available, w e imputed all fertility t o the married population.

A s with the multi-state life table program, t h e program MSLESLIE.PRC is the s a m e as f o r t h e standard Leslie matrix except t h a t the basic argument is the tran- sition matrix r a t h e r than t h e s c a l a r rate. The subdiagonal blocks are calculated from t h e multi-state life table using t h e expression hL, _*inv (hL,) and t h e f e r - tility blocks are similarly expressed using a matrix version of the standard formu- l a (see, f o r example, Keyfitz 1985, p. 206).

Projections were generated in five y e a r intervals f o r t h e 50 y e a r s following 1981. The projected v e c t o r is l a r g e (72 elements) and w e have taken as examples t h r e e summary measures: total population, number 6 5 + , and number 65+ and not married (Table 2). all f o r Canadian females. Declining numbers are due to o u r using fixed 1980-82 b i r t h rates t h a t are considerably below replacement. No account i s taken of immigration. Figure 1 shows t h e proportion of 65+ women who are not married when four alternative masks are imposed on t h e transition r a t e s .

Variances in the multi-state model

The multi-state analysis considered up t o t h i s point, like traditional life table analysis, i s based e n t i r e l y on deterministic t h e o r y . While t h i s a p p r o a c h may b e justified on a number of grounds, t h e r e are occasions when a stochastic analysis might b e considered more a p p r o p r i a t e . In p a r t i c u l a r , w e may b e i n t e r e s t e d in a stochastic model in cases where a ) t h e population size i s small or b) t h e transition intensities are thought to b e random.

I t may be helpful at t h i s point to d e s c r i b e briefly what i s meant by a stochas- t i c population model. Consider t h e matrix p, in t h e f i r s t section. I t specifies t h e probability t h a t a n individual at a g e

z

and state j will be in state i a f t e r one time period, given t h e underlying Markovian transition rates. When w e c o n s t r u c t t h e multi-state life t a b l e , however, w e treat p, as a matrix of proportions.

This assumption, t h a t probabilities may b e t r e a t e d as proportions, i s usually invoked in t h e construction of t h e single state life table, and i s justified by assum- ing t h a t t h e population i s sufficiently l a r g e s o t h a t random deviations from t h e ex- pected values, at l e a s t on a p e r c e n t a g e basis, may be ignored. If we use a s e r i e s of transition rates t o make predictions a b o u t t h e behavior of s m a l l populations sub- j e c t to c e r t a i n fixed r a t e s , however, w e should recognize t h a t population size and i t s distribution will b e random v a r i a b l e s with nonzero v a r i a n c e and i n t e r p r e t t h e p r o j e c t e d counts as e x p e c t e d values of t h e s e random variables. The idea of a st*

c h a s t i c population model i s to provide a convenient mechanism f o r gauging t h e size of t h e resulting variances. A g e n e r a l class of s t o c h a s t i c population models are t h o s e which may b e r e p r e s e n t e d by Galton-Watson branching p r o c e s s . The most im- p o r t a n t a s p e c t of t h i s class of models i s t h e assumption of independence among in- dividuals.

Pollard (1966) h a s provided a convenient method of computing v a r i a n c e s and c o v a r i a n c e s of a g e n e r a l Galton-Watson branching p r o c e s s . Of p a r t i c u l a r i n t e r e s t i s t h e f a c t t h a t a n extended Leslie matrix t h a t p r o j e c t s second and h i g h e r moments of t h e population a g e distribution may b e constructed using t h e d i r e c t or K r o n e k a r product. This r e s u l t also e x t e n d s t o t h e multi-state Leslie matrix, even if t h e rates themselves are random (Pollard 1983). I t i s t h u s possible to c o n s t r u c t a r a t h e r simple program, again using GAUSS, t h a t will provide estimates of v a r i a n c e f o r a l a r g e class of multi-state projections.

The method developed by Pollard r e q u i r e s t h e construction of f o u r matrices, T,

F,

8 , and

M,

whose matrix p r o d u c t i s t h e required extended projection matrix.

The p r o c e d u r e MSVARI.PRC c a l c u l a t e s this extended projection matrix given t h e

same transition intensity matrix discussed above. While a n initial version of t h e p r o c e d u r e actually c r e a t e d t h e s e f o u r matrices and t h e n multiplied them, i t t u r n e d out to b e more convenient to partition t h e f o u r matrices and c a r r y o u t t h e multi- plication p a r t i t i o n by partition. The r e a d e r should h a v e l i t t l e t r o u b l e following t h e code a p p e a r i n g in Appendix 11; however. i t may b e helpful to briefly discuss what e a c h matrix i s designed to do.

Matrices F and M t r a n s f o r m moments about t h e origin to falling f a c t o r i a l mo- ments and back again. They are lower t r i a n g u l a r , and at least f o r t h e t w o moment c a s e , e a s y to c o n s t r u c t . The matrix B transforms f a c t o r i a l moments of t h e set of in- itial states into f a c t o r i a l moments of a set of intermediary random variables, where e a c h intermediate random v a r i a b l e c o r r e s p o n d s to one t y p e of transition t h a t will b e o b s e r v e d in a given period. The transition to d e a t h and t h e identity transition must both b e included explicitly so t h a t t h e sum of transition probabilities associ- a t e d with e a c h state i s one. If P i s t h e block diagonal matrix with t h e i t h block be- ing t h e v e c t o r of transition probabilities associated with t h e i t h s t a t e , then w e c a n write B simply as a block diagonal matrix with blocks P and P.*.P

(.*.

^{i s used} to denote t h e K r o n e k a r or d i r e c t product).

The number of individuals in e a c h state at time t +1 will b e a l i n e a r t r a n s f o r - mation, denoted by t h e matrix Q, of t h e intermediate random v a r i a b l e s which in t u r n are derived from t h e number of individuals in e a c h state at time t . If only transitions among k explicit states are possible t h e n t h e r e will b e .up to k possible transitions f o r e a c h state (including t h e identity transition). The number of i n t e r - mediate random v a r i a b l e s , a n d t h u s t h e number of columns of

9,

will b e n o l a r g e r than k*k. In t h e case where a l l transitions are included t h e Q matrix will consist of k k-dimensional identity matrices concatenated horizontally. If one state (e.g.

d e a t h ) i s included implicitly, t h e n e a c h p a i r of identity matrices .will b e s e p a r a t e d by a column of z e r o s . In t h e Leslie projection model, as a r e s u l t of t h e p r e s e n c e of b i r t h p r o c e s s e s , Q will b e more complicated. In e i t h e r c a s e , t h e matrix T relating t h e f i r s t t w o moments of t h e intermediate and final random v a r i a b l e s will b e block diagonal with blocks Q a n d Q.*.Q.

Since all of t h e above t r a n s f o r m s are l i n e a r and s i n c e t h e matrices TI

M I

and F are fixed by t h e set of possible transitions, generalizations to changing and even independent random t r a n s i t i o n rates i s straightforward. In t h e iid random c a s e , f o r example, w e simply r e p l a c e B with i t s expectation, which r e q u i r e s t h a t w e specify E (P. *.P) a n d E (P).

The p r o c e d u r e MSVARI c o n s t r u c t s t h e r e q u i r e d matrix product given t h e a r - guments

P ,

Q , and optionally

E(P.*.P).

Table 3 provides a simple example based on a Leslie projection matrix. The p r o c e d u r e MSPQT g e n e r a t e s t h e requisite P and Q values from a d i s c r e t e time Markov transition matrix. An applications program then uses t h e s e p r o c e d u r e t o c o n s t r u c t a multi-state life t a b l e with variances from t h e Canadian marital s t a t u s d a t a a p p e a r s with t h e o t h e r p r o c e d u r e s in Appendix 111; a few e x t r a c t e d r e s u l t s may be found in Table 4.

In t h e simple c a s e with a fixed set of age specific rates and a fixed initial po- pulation t h e r e i s no need t o r e l y on t h e above a p p r o a c h because t h e variance- covariance matrix at a g e

z

may b e calculated directly from t h e 1, matrix using t h e multinomial distribution. The branching p r o c e s s i s needed f o r a model t h a t allows t h e r a t e s t o be stochastic o r includes renewal by b i r t h . W e are unable to estimate v a r i a n c e s using t h e (18 X 4) block-Leslie projection matrix described in t h e preceding section because t h e projection matrix would b e 72+72 X 72

=

5256 X 5256, which i s t o o l a r g e t o b e handled efficiently, hence t h e smaller ma- t r i x of o u r example.

The assumption t h a t individuals are independently s u b j e c t t o t h e observed transition probabilities r e s u l t s in v a r i a n c e s t h a t are too small in application t o l a r g e populations. The Pollard method initiates but by no means exhausts t h e study of v a r i a n c e in a branching p r o c e s s .

Figure 1

Projected proportion of 65+ Canadian women who a r e not married (after abstracting from specified transitions).

Base Year (1 981=0)

t No Rerna 0 No Div _A No Wid

Table 1

Extract from life table by a g e and marital status.

Canadian Females, 1981

Table l a :

No

Mask

Tran. Originating Status Orlginatlng Status

Age Status S M W D M W D

Table l b :

No

Remarriage Mask

Tran. Origlnatlng Status Orlglnatlng Status

Age Status S M W D M W D

Table 2

Projected Female Population (1000s) Base Run

Year Total 65

+

6 5 + and Not Married

T a b l e 3

Illustration of Stochastic Leslie Matrix

T a b l e 3a

Intermediate transition probability matrix

Transition Type

Age group

1 2 3

T a b l e 3 b

Intermediate transition counts matrix

Age Transition Type Group 1 2 3 4 5 6

T a b l e 9c

Stochastic Leslie Projection Matrix

Note: The original 3 X 3 Leslie matrix a p p e a r s in t h e upper left c o r n e r of the stochastic matrix.

Table 4

Expectations and variances of L 65 from a condensed stochastic life table by age and marital status.

Canadian Females, 1981. Radix age 20.

1,

Trans. Originating status

Status Moment Single Married Widowed Divorced

Single M

v

Married M

v

Widowed M

v

Divorced M

v

REFERENCES

Edlefson, L.E. and Jones, S.D. (1986) G A m . Copyright software.

Feeney, G . (1970) Stable a g e by region distributions. Demography 7:351-348.

Keyfitz, N. (1985) Applied Mathematical Demography. N e w York: Springer- Verlag

.

Pollard,

J.H.

(1966) On t h e use of t h e d i r e c t matrix product in analyzing c e r t a i n stochastic population models. Biometrika 53(3):397-414.

Pollard, J .H. (1973) Mathematical Models for t h e Growth of Human P o p u l a t i o n s . Cambridge, England: Cambridge University P r e s s .

Rogers, A. (1975) I n t r o d u c t i o n t o Multiregional Mathematical Demography. N e w York: John Wiley & Sons.

Schoen, R. (1975) Constructing increment-decrement life tables. Demography 12:313-324.

Appendix I

@

Program: MSLIFE. PRC

@

@

This program produces a procedure which constructs a multi-state life table from the argument "vcm", a vertical concatenation of age specific transition rates. It

is

assumed that there are

5

year age groups. Transitions occuring before the radix age are ignored.

@

proc mslife(vcm,radage);

local k, n, m, ages, ix, id, vcp, id, vcl, vc11, vc112, vct, vce,

i ,

ria, rill;

k=cols (vcm)

; @

k

is

number of states

@

vcm=trim(vcm, k*radage/5,0)

; @

Discard data below radix age

@

n=rows (vcm)

;

m=n/k;

@

m

is

number of age groups

@

ages=seqa (radage, 5 , m)

; @

Creates vector of ages

@

ages=vec

(

(ages. *ones (1, k)

)

'

⁾ ; @

Ages vec matched to VCM rows

@

ix=reshape (seqa ( l , l , n) , ^{m, k)} '

; @

ix

is

matrix used for indexing

@

@

VC matrices

@

id=eye (k)

; @

id

is

k dim identity matrix

Cd

vcp=vcm; vc l=vcm; vc 1 l=vcm;

@

Matrices set to correct size

@

VCE=VCM; VCT=VCM;

clear rill;

@

Loop over ages to create px, lx, 5Lx

@ i=0;

do until i==m; i=i+l;

riB=ixC.

, i l ;

vcpCri0,. l=inv(id+2.5*vcmCri0,. I >*(id-2.5*vcmCri0, .

^{1 > ;}

if i==l;

vclC ria, . ^{I =id;}

else;

vclCri0,. l=vcpCrill, . I*vclCri11,. I

;

vcllCril1,.

I

=(vclCri0, . ^I +vclCrill,.

1 )*5/2;

if i==m;

vcllC ria, . ^{I =vclC ria,} . ^I

^*5/2;

endif;

endif;

rill=ri0;

endo;

@

Loop over ages to create Tx and Ex

@

vc112=reshape (vcll, m, k*k>

; i=0;

do until i==m; i=i+l;

ri@=ixC. , ¹¹

;

vctCri0,.

l=reshape(sumc(trim(vcll2,

i-I,@)), k , k);

vceC ria, . ^{I =vctC ria,. I} *inv(vclC ria, . ^I

^{) ;}

endo;

@

Return altered matrix, close the procedure, and save to disk@

r e t p ( a g e s " v c p " v c l " v c 1 1 ~ v c t ~ v c e ~ ;

endp;

save mslife;

@

Program: MSMASK. PRC

@

@

This program produces the procedure MSMASK and a few masks of

interest for the analysis of marital status which allow one to zero certain transition probabilities according to a k by k mask of

ones and zeros; ones indicate that a particular element should NOT be altered. A zero along the diagonal will cause the death rate in that group to be set to zero. In all cases the diagonal elements will be adjusted so that the sums along the columns are equal to the death rates.

@

proc msmask (vcm, mask

1 ;

local n, k , m, ix, nmask, ria, i, mx, dthx, mx;

n=rows (vcm)

;

k=cols (vcm)

;

m=n/k;

ix=reshape (seqa (l,l, n) , m, k) '

;

@

Apply mask and recalculate diagonal elements

@

i=0;

do until i==m; i=i+l;

ri@=ixC., il

;

mx=vcmC ria, .

^{1 ;}

dthx=sumc (mx). fdiag(mask1;

nmask=mask. f (-eye <k) +1)

;

mx=nmask. Xmx;

vcmC ria,

. I

==+eye (k) . ^X (dthx-sumc (mx)

1 ;

endo;

@

Return altered matrix, close the procedure, and save to disk@

retp (vcm)

;

endp;

@

Each mask sets certain transition probabilities to

zero by applying a 4 by 4 "mask" to the data. Here we provide four possible masks which may be applied to marital status data ordered as: single, married, widowed, divorced. The base mask will have no effect.

@

Base=ones (4,4

1 ;

let NoReMaC4,41=1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 ; let NoDivC4,41 =1 1 1 1

1 1 1 1 1 1 1 1 1 0 1 1 ; let NoWidC4,41 =1 1 1 1

1 1 1 1 1 0 1 1 1 1 1 1 ; let NoDeaC4,41

=0

1 1 1

1 0 1 1 1 1 0 1 1 1 1 0 ;

save msmask, NoDea, NoWid, NoReMa* NoDiv* base;

@

Program: MARRIAGE. CAN

@

@

Multi-state life tables can be created using marital status data from Canada. Elements of transition matrices may be set to zero using a masking procedure. Choose an appropriate data set by changing the two lines following this comment.

@

loadm cf8;

vcm=cf

8 ;

radage=20

;

k=4

;

n=rows (vcm)

;

m=n/k;

loadp mslife, msmask;

loadm base, norema, nodiv, nowid, nodea;

@

Certain transition probabilities may be set to zero by applying a 4 by 4 "mask" to the data. Here we provide five possible

masks: all rates (base) no remarriage (NoReMa), no divorce (NoDiv), no widowhood (Nowid), and no death (NoDea). To select one simply use the appropriate condition name in the call to

MSMASK.

Combination masks may be specified simply by using

.AND (e.g. NODEA .AND ITOWID to omit death of person and of spouse

)

.

@

vcm=msmask (vcm, base

) ;

@

Call multi-state life table procedure

@

lft

=

mslife (vcm, radage)

;

ages=lftC . , 11

;

vcp=lftC . ,2: 51

;

vcl=lftC .

^,6:

⁹¹

^;

vcll=lftC . ^{, 10: 131}

^;

vct=lftC . , 14: 171

;

vce=lftC . ,18: 211

;

@

Control printing

@

format 10,3;

fn rd (x) =O. 001*round (1000*x)

;

keepages=20

1

65;

i=0;

do until i==rows (ages)

;

i=i+l;

if sumc (ages[ i, 11 . ==keepages) >O;

rd(agesCi,.l"vclCi,.I"vceCi,.I~;

if i/4==trunc (i/4)

;

print; endif;

endif;

endo;

Appendix I 1

@

Program: MSLESLIE. PRC

@

@

Produces procedure that returns Leslie matrix for the multi- state projection of female population with fertility rates cvasfr and transition rates cvm. Age structures of rates are those for Canada: population data has 16 5-year age transition data has 17 5-year age groups; Fertility

is

age specific

fertility 10-49 in single year age groups. It

is

assumed that these rates apply to married people, and that there

i s

no childbearing outside of marriage.

@

loadp mslife;

proc msleslie (cvasfr, cvm)

;

local n,m, k , x, vc11, ix,asfr5, vmf, les, i,ri0,ri11;

n=72

;

m= 18

;

k=4

;

x=mslife(cvm,

0);

vcll=xC . , ^{10: 131}

^;

ix=reshape (seqa (l,l, n) , m, k) '

;

asf r5=zeros (mf 5 , 1 >

;

asfr5C 10: 49,ll =cvasfr;

asf r5=sumc (reshape (asf r5 , m, 5) '

) ;

vmf=zeros (m, kfk)

;

vmfC . ^{,21 =asfr5;}

vmf=reshape (vmf, mfk, k)

;

les=zeros (n, n)

; i=0;

do until i==m; i=i+l;

ri0=ixC. ,

il ;

if i>1;

@

Mortality projections using ratio of 5Lx for consecutive years

@

@

Fertility projection using product of infant survival (5L0/5>, average 5 year fertility for member of xth age group (5Fx+

5Fx+5*5Px>/2, and fraction of children female (1/2>

@

endif;

rill=ri0;

endo;

retp (les)

;

endp;

@

Program: MSPROJ. CAB

@

@

This program uses the MSLESLIE procedure to project the female population by age and marital status in Canada for 50 years in 5 year steps. Data

is

aggregated in a variety of ways and numbers are plotted. For notes on masking see marriage.can and msmask.prc.

@

loadm ccf8,cf8,casfr;

#include "msleslie. prc"

;

loadp msmask;

loadm base, norema, nowid, nodiv, nodea;

cvm=msmask (cf

8,

base

1 ;

popmat=zeros(72,11>;

popv=vec (ccf 8'

) ;

les=msleslie (casfr, cvm)

; i=0;

do until i==11; i=i+l;

popmatc . ,

il

=popv;

popv=lestpopv;

endo;

format 8'2;

@

Create summary measures of population projection, p-t

:

total female population by year

p-o

:

number of females older than 65

p-o-nm

:

number of females older than 65 and not married

@

years=seqa (0,5,11>

;

ages=vec (seqa (0,5,18). tones(1,4)

) ;

states=reshape (seqa (1,1,4), 72'1)

;

old=ages. >=65;

notmar=states. /=2;

all=ones (72'1)

;

p ~ o ~ n m = s u r n c ( p a c k r ( m i s s ~ 0 1 d

.and notmar,0)"popmat>>;

p-o=sumc(packr(miss(old,

0)"popmat));

p-t=sumc (packr (miss(al1, O)Mp~pmat)

) ;

print

years"trirn(p-tMp-o"p-0-nm,

1'0);

Appendix

_{I 1 1}

@

Program: MSVARI . ^PRC

^@

@

This program creates a procedure MSVARI which returns the projection matrix for the first two moments of a multi-type branching process. Three arguments are used.

P

is

the block diagonal combination of the column vectors

p(i> of probabilities of possible transition sets from the ith state.

Q

is

the matrix of counts for each transition state. The ith row has the number of counts that should be added to the

ith state for each transition type.

PKP is the expectation of the own direct product of p. Set to zero unless transition rates are themselves random so that PKP/=P. *. ^P

^@

proc msvari(P,Q,PKP>;

local sts, trn, i, T11, B11, T22, B22, F 2 1 , All, A12, A22, A21, A;

sts=rows

(q) ;

trn=cols(q>

;

Tll=Q; T22=Q. X. Q;

B11=P; B22=P. X. P;

F21=zeros (stsXsts, sts)

;

i=8

;

do until i==sts; i=i+l;

F211: (1-l>fsts+i, il=-1;

endo;

M21=zeros (trnXtrn, trn)

;

I=@; do until i==trn; i=i+l; m21C (i-l>ftrn+i, il=l; endo;

All=Tll*Bll;

A21=T22fM2lfBll+T22tB22tF21;

A22=T22XB22

;

A12=zeros (sts, stsfsts)

;

A=All"A12: A2lWA22;

retp (A)

;

endp;

save msvari;

@

Program: MSPQT. PRC

@

@

This program creates a procedure MSPQT which takes as an

argument a Markov transition intensity matrix and produces the P and Q (transposed) matrices that are needed by the MSVARI

procedure.

@

proc mspqt (mu)

;

local k, id, p, p2, bigp, bigqt,

i ,

ik, vp, pak;

k=rows (mu)

;

id=eye (k>

;

p=inv(id+2,5tmu> t (id-2.5tmu)

;

p2=p: (-sumc (p>

'

+1>

;

bigp=zeros(kt (k+l), k >

;

bigqt=bigp;

i=0;

do until i==k; i=i+l;

ik=seqa( (k+l>t(i-l)+l, 1, k+l);

bigpC ik, 11 =p2C. ,

^{il ;}

bigqtC ik, . ^{I =eye (k) : zeros(1, k);}

endo;

vp=miss(vec(p2>,

0) ;

pak=packr(vpWbigp-bigqt);

pak=trim(pak' , l,0) '

;

retp (pak)

;

endp;

save mspqt;

@Program: LESLIE. EXP @

@ T h i s program i l l u s t r a t e s t h e P o l l a r d method f o r c a l c u l a t i n g

v a r i a n c e s i n a s i m p l e 3x3 L e s l i e m a t r i x a p p l i c a t i o n . The o r i g i n a l L e s l i e m a t r i x may be r e c o n s t r u c t e d by t h e p r o d u c t q1p.

(3

l o a d p m s v a r i ;

l e t pC6,31= 1 0 0

0 0 . 1667 0

0 0 . 1667 0

0 0 . 3 3 3 3 0

0 0 0.6667

0 0 0 . 3 3 3 3 ;

l e t q C 3 , 6 1 = 0 1 1 0 1 2

1 0 0 0 0 0

0 1 0 1 0 0 ;

format 5 , 2 ;

f n r d ( X I = @ . 011round ( 1 0 0 1 ~ ) ; r d c m s v a r i ( p , q , 0 ) 1 ;

rQ

Program: MARVAR. CAN

@

@

Multi-state life tables with variances are created using

Canadian marital status data. Elements of transition matrices may be set to zero using a masking procedure. For notes

on data and masking see marriageecan

@

loadm cf8;

vcm=cf

8 ;

radage=20

;

k=4

;

n=rows (vcm)

;

m=n/k;

#include "mslifeva.prc";

loadp msmask;

loadm base, norema, nowid, nodiv, nodea;

vcm=msmask (vcm, base)

;

@

Call multi-state life table procedure

@

lft

=

mslifeva(vcm,radage>;

ages=lftC . ^{, 11}

^;

states=lftC . ^{, 2}

^{31 ;}

vcl=lftC . ^{,4: 71}

^;

@

Control printing

@

format 10,3;

fn rd

(x)

=O. 00ltround (1000tx)

;

keepages=20 : 65

; i=0;

do until i==rows (ages)

;

i=i+l;

if s u m (ages[

i ,

11 . ==keepages1

>0;

rd(agesC

i ,

. ^I "states[

i ,

. ^{I -vclC}

^{i ,}

.

^{I ) ;}

if i/20==trunc(i/20>; print; endif;

endif;

endo;

Appendix I V

@Program: MAKEDATA. CAN@

@Makes the data needed for multi-state applications programs and saves to disk in the form of GAUSS matrices. Data are for

Canadian females, 1981. Elements of cf8 are transition

intensities by age

(0-89

in 5 year groups) and marital status (Single, Married, Widowed, Divorced). Elements of ccf8 are the 1981 population by age and marital status for the same age and marital status groups.Elements of casfr are single year age specific fertility rates from ages 10-49 for all women.

@

let cf8[ 72,41=

0.0021

0 0 0

0 0 0 0

0 0 0 0

0 0 0 0

let casfrC 40,1l=

0 0 0

0.0002 0.0012 0.0045 0.0122 0.0234 0.0360 0.0494 0.0637 0.0805 0.0970 0.1110 0.1216 0. 1310

let ccf8C 18,41=

save ccf8,cf8,casfr;

Appendix V

The r e s u l t s f o r t h e b a s e r u n of e a c h p r o g r a m may be r e p r o d u c e d i n t h e f o l l o w i n g way. F i r s t , e n t e r e a c h of t h e p r o g r a m s i n ASCII f o r m a t i n c l u d i n g t h e d a t a c r e a t i o n p r o g r a m a n d t h e f o l l o w i n g s h o r t p r o g r a m which may be u s e d t o c r e a t e t h e p r e c o m p i l e d p r o c e d u r e s :

@ Program: MSSTART. CAN @

@ T h i s p r o g r a m w i l l create t h e r e q u i r e d d a t a matrices a n d p r e c o m p i l e d p r o c e d u r e s a n d s a v e t h e m t o d i s k .

@

# i n c l u d e " m a k e d a t a . c a n " ; ;

# i n c l u d e " m s l i f e . p r c " ; ;

# i n c l u d e "msmask. p r c " ; ;

# i n c l u d e " m s v a r i . p r c " ; ;

# i n c l u d e " m s p q t . p r c " ; ;

Each p r o g r a m s h o u l d b e i n a s e p a r a t e f i l e , a n d t h e f i l e name s h o u l d be t a k e n f r o m t h e f i r s t l i n e o f t h e p r o g r a m . Then, f r o m DOS t y p e :

GAUSS <RTN>

RUN MSSTART.CAN < F 4 > < F 2 >

You may t h e n r u n e a c h o f t h e f o u r a p p l i c a t i o n s p r o g r a m s by t y p i n g : RUN MARRIAGE.CAN < F 4 > < F 2 >

R U N MSPROJ. CAN < F 4 > < F 2 >

R U N LESLIE.EXP < F 4 > < F 2 >

RUN MARVAR.CAN < F 4 > < F 2 >