- Email: [email protected]
Simulation of micropopulations in epidemiology: Tutorial 1. Simulation: an introduction
Simulation of micropopulations in epidemiology: Tutorial 1. Simulation: an introduction A series of tutorials illustrated by coronary heart disease models Eugene Ackerman* National Micropopulation Simulation Resource, Health Computer Sciences, Box 511 UMHC, University of Minnesota, 420 Delaware Street SE, Minneapolis, MN 55455, USA
(Received 31 August 1993;accepted 8 September 1993)
Abstract
This is the first in a series of tutorials concerning the simulation of micropopulation models for use in epidemiological research. The series emphasizes the techniques that are in use at the National Micropopulation Simulation Resource at the University of Minnesota. This initial series of tutorials uses chronic disease models, specifically coronary heart disease models, to illustrate the principles and methodologies employed. The series is intended for persons having some familiarity with electronic, digital computation and desiring to know more about a specific class of simulations, namely those that apply Monte Carlo techniques in the simulation of micropopulation models as a part of biomedical research. It is the author’s hope that these tutorials will be particularly useful for students in medical informatics; they could also be of interest ,to professionals in both the computational and the biomedical sciences. Keywords: Simulation;
Micropopulation;
Risks; Stochastic model; Chronic disease model; Monte Carlo
1. Simulation Simulation is an extension of the concept of modeling complicated systems by simpler ones; this has been an important technique in the natural sciences for many years. For example, physicists and engineers work with mathematical representations of frictionless systems to learn about the behavior of real objects. Biologists work with animal models to represent human organ systems and * Corresponding author. E-mail: [email protected]. umn.edu.
with mechanical models to represent circulatory dynamics. The use of models in no way obviates the need for careful experimental observations; in fact, quite the opposite is true. Nonetheless, the models, by their very simplifications, make it easier for humans to understand natural phenomena. For many years, there was a widespread suspicion of mathematical models within biological communities. This distrust was frequently based on the belief that the simplifications necessary to arrive at a mathematically analysable model
0020-7101/94/~07.00 0 1994 Elsevier Science Ireland Ltd. All rights reserved SSDI 0020-7101(93)00939-F
230
E. Ackerman / Int. J. Biomed. Comput. 36 (1994) 229-238
overlooked too many key features of the real system. All too often that was a correct interpretation. However, since the advent of computer-based analyses and simulations, over-simplifications are no longer a necessary part of modeling. Computer-based simulation has gained respectability within the natural sciences and the social sciences. The applications of simulation techniques are so widespread, and the models so diverse, that it would not be appropriate to catalogue them in this tutorial. Neither would it be very instructive, since one could read a list of applications with no understanding of the simulations involved. Rather, greater detail is presented here about a limited group of simulation models. In some discussions of simulation, it is convenient to distinguish between deterministic models and stochastic ones. In the former, there is no uncertainty concerning the outcomes of the simulations, whereas in the latter the outcomes are distributions of events and of uncertainties about those distributions. In actual practice, the ideal of a purely deterministic or a purely stochastic model is never met; all actual models have some elements of both. This set of tutorials is concerned with models that are primarily stochastic in nature. The models presented are called micropopulation models in that each individual in the population is followed through a sequence of time periods. The models discussed are ones that are represented by difference equations rather than differential ones; the outcomes of these models are determined using Monte Carlo techniques. These terms are defined and examined more fully in this first tutorial. Although the restrictions in the preceding paragraph exclude most uses of simulation in biomedical research, the possible areas of the models described remains very broad. For pedagogic purposes, it seems suitable to limit such applications still further. Accordingly, this set of six tutorials describes the use of simulation in epidemiological studies of coronary heart disease. This first tutorial is meant primarily as an introduction to this type of simulation study and the last of the six tries to place it in a broader context. The four tutorials in between consider various aspects of these simulations. A frequent use of modeling in epidemiology is
concerned with macropopulation models that are primarily statistical in nature. These will not be discussed further here. The interested reader can find a comprehensive bibliography with comments in Ref. 1. The studies described in the current set of tutorials are on-going, and it is not possible to predict their future directions. Although emphasis is placed on completed studies, none the less studies in progress are mentioned. The results of the continuing studies should appear eventually in journal article format. Additional tutorials, however, may focus on other applications of the techniques employed at the National Micropopulation Simulation Resource. 2. Chronic disease -
coronary heart disease
Epidemiology originally meant a study of epidemic occurrences of disease. Epidemic, in turn, meant that it was outside of, or in excess of, the number of cases of that disease expected in the given population. Epidemics might be due to infectious agents or to transitory changes in the environment. Over the years the term epidemiology has been changed and extended to mean the study of the occurrence and etiology of disease in a population. The present series of tutorials is concerned with simulations of models of interest in epidemiology. In particular, this series has been restricted to models of chronic diseases, that is to ones that continue over an extended period of time. The diseases considered are non-infectious ones. For simplicity, only non-infectious diseases associated with the heart and the coronary circulation are considered. These are described by the label coronary heart disease and will be referred to by the letters CHD; the abbreviation is used in this fashion in numerous references, one of which is Ref. 2. A review of coronary heart disease epidemiology can be found in Ref. 3. For these simulation studies, it is convenient to consider a population as composed in part of persons who have not shown symptoms of CHD. Others in the population may have heartassociated pains called angina with or without any other evidence of CHD. The population will also be considered to include individuals who have
E. Ackerman /ht.
J. Biomed. Comput. 36 (1994)
experienced one or more non-fatal myocardial infarcts. The latter are occurrences in a region of the heart muscle myocardium in which the local circulation fails resulting in a permanent loss of function in that region of the myocardium. Where convenient, the abbreviation MI will be used to designate myocardial infarct. Other types of CHD are well known, but will not be considered in this series of tutorials. Rather, it will be assumed that all individuals are free from CHD except for possible angina and previous MIS. Some MIS are called silent because the individual does not report any symptoms; silent MIS are known only from autopsy data. In contrast to these, all of the MIS considered in the models presented are so-called acute MIS (AMIs). The latter may be fatal or non-fatal and may or may not involve hospitalization. Amongst the other events that may happen to a simulated individual are ones that involve death. Death may occur due to causes other than CHD; this type of event is referred to by the letters OCD. Death due to CHD (CHDD) may be preceded by, or the consequence of, an AM1 or simply occur suddenly with other etiology. The latter is referred to as sudden CHD death. The probability of experiencing CHD is highly age-dependent. Age may be specified as a risk, or the population simulated may be restricted to individuals within a given age range. In addition, at least in the past, the probability of CHD has been found to be quite different for men than for women. Thus, it may be necessary to carry out the simulations for sub-populations of only one sex. Heterosexual populations can be resolved into two separate ones for simulation purposes, or the simulation model may use the individual’s gender as a pointer to a sub-model. 3. Goals for coronary disease models The selection of a particular type of simulation model and of a set of simulation techniques is highly dependent upon the goals of the investigator. One goal that is emphasized in this set of tutorials is predicting the changes that might be anticipated as a result of population-based intervention programs. Another goal might be the evaluation of various risk factors, both in their
229-238
231
quantitative importance and also in the relative strength of their effects. A quite different approach is to seek a model that allows one to incorporate knowledge of the physiological processes underlying the pathologies modeled. By and large, in the area of chronic disease epidemiology, such models have not been widely used. By contrast, in infectious disease epidemiology, physiologically-based models are more prevalent [4]. However, CHD models that incorporate even a small approach to representing physiological processes are more intellectually satisfying than ones that make no attempt to do so. Another goal, albeit one that is often not mentioned explicitly, is to have a model that allows an easy estimate of the prevalence or the incidence of CHD-related events. The ease may depend upon the availability and user-friendliness of packaged statistical programs. Closely related to that goal is the one of simplicity. If two models are equally successful in other regards, the simpler one appears preferable. (The criterion of simplicity is sometimes called Occam’s razor.) Other goals are also important. One concerns the design of new surveys and new clinical trials. The use of a predictive model can help select the items to be determined, the population size and the duration of the surveillance. Similarly, model simulations can suggest the most critical parameters to monitor. In turn, the new studies can help to test the selected model. A frequent need is the communication of the results of a study to others. The use of models and of micropopulation simulations can make this communication easier. Thus the model and the simulation become a conceptual framework rather than a separate investigational tool. Goals will be referenced repeatedly in the discussions of specific models. By and large, a given model may be optimal for one of the goals cited in the preceding. Few, if any, models are well suited to all of the aforementioned goals. 4. Discretized stochastic micropopulation
models
In CHD studies, emphasis is placed on the occurrence of one or more of several possible events. Events, such as the occurrence of an MI, are dis-
E. Ackerman
232
/ ht.
J. Biomed.
Crete ones. Individuals experiencing such an event are said to move to a new state. Population models frequently involve simulations in which the fraction of the total population in a given state is considered. A model of that nature is said to be a macropopulation one. For such models, no attempt is made to keep track of an individual population member or an individual’s past history. On the other hand, when sufficient heterogeneity is present in the population, it is easier to follow each simulated individual. This type of model is called a micropopulation one. Only micropopulation models are considered in detail in this series of tutorials. The states of a model of that nature are illustrated in Fig. 1, which shows a flow diagram for a stochastic, micropopulation model that could be simulated using suitable techniques. The model is simulated in general for a limited period of time, for example 7 years. The length of the simulation is referred to as an epoch. During a simulation, all individuals who can experience some event, for example an acute MI, are followed either until they experience an event or until the end of the epoch, which ever occurs first. By and large, there is a low probability of a given type of event, so the largest group of simulated individuals
I
Death
from
Coronary Heart
Disease
No Previous History of Myocardial Infarct
Fig. I. Flow diagram of a model of coronary heart disease. This is an extremely simplified model of three possible events that may occur to a healthy individual. This and other simplified forms are used in several of the following tutorials.
Comput.
36 (1994)
229-238
will start with no experience of CHD. and continue that way throughout the epoch. In micropopulation models, it is possible to represent the probabilities of events as continuous functions of time which may be described by differential equations. However, in the presence of population heterogeneity, it is easier to deal only with discrete events and discrete units of time. That is true for all of the models to be treated in this and following tutorials. If time is discretized, then there is no meaning to trying to represent the events as occurring within a part of the basic time interval. Even though the program representing a micropopulation model may consider the individuals serially (that is, one after the other), it is not possible for events occurring during that interval to influence changes occurring during that interval. Thus the models are well described by difference equations; values of interest are determined by values at the previous time interval. The difference equations may be represented explicitly or may be described algorithmically as part of the simulation program. 5. Monte Carlo and other simulation techniques Models of the nature discussed here are readily simulated by Monte Carlo techniques in which pseudo-random numbers are used to decide if a given individual experiences a possible event in a discretized time period. Fig. 2 illustrates the fashion in which a uniform random variate can be used to select amongst several possible discretized outcomes. Effectively, this method transforms a uniformly distributed random variate to one that is distributed in the fashion appropriate to the simulation being considered. (It is also possible to transform a uniformly distributed random variate to any arbitrary continuous distribution using a technique similar to that illustrated in Fig. 2. Selection of a level or value based on a continuous probability function is used in some micropopulation simulations.) For illustrative purposes, assume that there are two possibilities from which one is to be selected. It is also necessary to assume that it is possible to compute the probability of each possibility. (The sum of the two probabilities must be 1.0.) Then a
E. Ackerman/Int.
J. Biomed. Comput. 36 (1994)
Outcome
I L-u t Event A
No
Change
0.0
0.2
Oi
Event B
0.6
Event C
08
1.0
Value of Random Variate
Fig. 2. Selecting outcomes using Monte Carlo techniques. The figure illustrates the process of choosing of one of four discrete options. First, a value is selected for a pseudo-random variate distributed uniformly in the range [O,l). If the value of the pseudo-random variate were 0. IS, i.e. below 0.4, the option No Change would have been selected, whilst if the value were above 0.8 the option Event C would have been chosen. Effectively this has transformed the uniformly distributed variate to a quadrunomial variate with probabilities 4 = 0.4, Pa = 0. I, P, = 0.3, and PC = 0.2.
random number uniformly distributed between 0.0 and 1.0 is selected. If the random number is smaller than the probability of the first possible outcome, that outcome is selected. Otherwise the second outcome is chosen. For more details please see Ref. 5. Where convenient, a wide variety of so-called packaged programs can be used to assist in the necessary statistical computations. Sets of packages that are referenced in successive tutorials include SAS [6], BMDP [7], and IMSL [S]. It is assumed that these are readily available to any group simulating stochastic micropopulation models representing the incidence and/or prevalence of CHD events in a given population or subset of the population. It is also assumed that any individual reading about such populations has at least an awareness of these specific packaged programs. The resulting model may be thought of as having at least three parts. One is comprised of the states and their interconnectivity as illustrated in Fig. I. This first part is sometimes referred to as a compartmental model since the states may be regarded as compartments of the population. However, the flow of individuals from one compartment to another is quite different than that of the usual meaning of compartmental analysis.
229-238
233
The second part of the CHD model is an algorithmic description of the procedures by which the probabilities implicit in Fig. 1 are determined. In some usages an algorithmic part is referred to as the statistical model, although this may be a quite misleading designator. A variety of algorithms are discussed in the next tutorial. These algorithms are described using one or another of the terms defined in Table 1; the notation used is modified from that found in Ref. 9. The third part of the CHD model of coronary heart disease presented is concerned with specific risk factors and their interactions. The specific risk factors included are chosen based on biomedical and epidemiological knowledge. The topics concerned with risk factors, their descriptions and their use in simulations are commented upon below in this tutorial and are developed more fully in following tutorials. 6. Field observations -
polychotomicity
The coronary disease models, since they are population-based and not driven by either known physiology or anatomy, must be verified by comparison with data on one or more real populations. The data are usually acquired by some type of survey or clinical trial. It is customary to refer to such observations as field data. Because of the fashion in which it is gathered, the data have far greater uncertainties than those associated with a laboratory experiment or with a controlled environment. The events of interest in coronary disease models occur relatively rarely. Thus the observations must be extended over a considerable period of time. Partly as a result of this, an individual is at risk for more than one possible type of event, several of which may exclude the others. This makes it harder to estimate the probabilities associated with each kind of event. Situations where only two outcomes are possible are called dichotomous, whilst those with multiple types of outcomes are designated polychotomous. One difficulty with all polychotomous models is estimating the parameters of the model. This is exacerbated by the mutual exclusivity of the events considered. One method of avoiding this limitation is to consider several subsets of the popula-
E. Ackerman I ht. J. Biomed. Comput. 36 (1994) 229-238
234
Table I Descriptors for dichotomous probabilities Name
Symbol
Event Density
Equivalent names
Definitions
Other relationships
event distribution pdfa death density
Pdr =f(t)dr Pdr is event probability during (1,~+ dt)
f(f) = dF(t)ldr f(‘,‘= -dS(r)/dt f(r)df s
s
= I
0
T
Cumulative Event Distribution
Survival Function
F(T)
S(T)
cdf” cumulative probability
survivorship death function
F(T) =
f(+ 0
-f(r)dr
S(T) =
s
T
F(T) = I - S(T) F(O)=0 F(a)
*
1
S(Q=l-F(T) S(0)= I S(m) * 0
Hazard Rateb
X(r)
force of mortality failure rate instantaneous death rate
f(l) = W)W) X(r) = -(dS(r)ldr]/S(r)
s
k(t) = -d[log(S(r))]ldr MO) =f(O)
T
Cumulative Hazard
A(T)
A(T) =
X(T)dr
0
A(r) = -logIS( T)l S(r) = e-*(n A(0) = 0 A(a) * m
“Abbreviations: pdf, probability density function; cdf, cumulative distribution function. bFor some purposes it may be desirable to use a slightly different definition for the hazard rate, namely A(T + 7,T) =/(T + 7)/S(T) which is the probability density of dying at time T + 7 conditional on survival to time T. This form is not used in these tutorials.
tion separately. Each smaller population contains the individuals experiencing no events or those individuals who experience a particular kind of event. To simulate the CHD model including the possibilities of different types of interventions, it is necessary to know not the probability of an event within the smaller population, but the probabilities within the entire population. Attempts to estimate these often fail for purely numeric reasons. Accordingly, an approximation is needed. The method described in the following section is valid only for events with a quite low probability. It permits estimations of the necessary parameters for polychotomous models, based on the more readily computed estimates from a group of dichotomous models. There are a number of different statistical terms that are used to describe the probability of an event. These are presented here for an individual in a dichotomous situation where the possibilities
are either an event that removes the individual from the population or no such event. These functions and distributions include probability density functions and cumulative distribution functions. The nature of these terms is such that if one is given, all of the others can be found, at least numerically. Table 1 presents a summary of these different functions; by and large, for the simulations of interest here, the cumulative distribution function (F) and the survivorship function (s) are the most important terms. A variety of expressions and rules can be used to describe the probability of the occurrence of an event in a given time interval, jI They are not commented upon further in this first tutorial, but are used throughout the series. 7. Multivariate risks factors
normalization of risk
It is sometimes pleasing to humans to assign a single cause to events that it is desirable to avoid.
E. Ackerman
/ht.
J. Biomed.
The cause is known as a risk factor. For coronary heart disease, it is necessary to include several different risk factors to simulate the events in the population. By and large, the risk factors have different units and are dimensionally incompatible. In order to combine them in mathematical expressions, it would be possible to multiply each risk factor by coefficients having the reciprocal dimensions. However, it then would be impossible to compare the relative importance of the risk factors since each coefficient would have different dimensions. In order to avoid that problem, one may use normalized dimensionless variates. A popular method of achieving this is to divide each risk factor by the population standard deviation. In order to further normalize the risk factors, it is also customary to subtract the population mean before dividing by the standard deviation. Thus a dimensionless, so-called normalized, risk factor represents the number of standard deviations that the actual risk factor is away from the population mean. This conversion may be expressed symbolically, letting yii represent the value of thejth risk factor for the ith individual. Then the corresponding dimensionless risk factor, x0 would be xv =
Yij - Yj OY,
where yj is the population mean for the jth risk factor and au. is its standard deviation. Although the unnormalized risk factors such as diastolic blood pressure and number of cigarettes smoked per day are always positive, the normalized values can be negative if the individual’s values are below the population mean. Dividing by some number such as the standard deviation is necessary to obtain dimensionless risk factors. However, the reference values need not be the population means; indeed, in some cases it makes more sense physiologically or epidemiologically if they are defined otherwise. The method of combining the risk factors to obtain an overall risk for a particular type of event such as an AMI is not dictated by any physiological or anatomic rule. Nonetheless, in the absence of any other indicator, it is customary for most models to use the exponential function of a linear combination of risk factors. This is a quite
Comput.
36 (1994)
229-238
235
arbitrary decision that is followed in most population studies. Provided the states and their interconnectivity, the algorithmic definition of the probabilities, and the form of the combination of the risk factors have been selected, it is usually possible to estimate the risk coefficients from the field data. As noted, this is most readily done by considering subpopulations with dichotomous outcomes. In the next few paragraphs, a method is described to estimate the risk coefficients for the entire population from the ones calculated for the subpopulations. In this tutorial, the normalized risk factors are represented by the symbol xii where the subscript i indicates an ordinal used to enumerate the individuals in the population and the subscript j, an ordinal used to enumerate the risk factors. Where appropriate the subscript k is used to indicate the particular outcome; the value k = 0 indicates no event. Since it is never possible to include all possible risk factors, the value j = 0 is used to indicate the combined influence of the unspecified risk factors. The value of xi0 = 1 is chosen arbitrarily for convenience in notation; it merely indicates a lack of knowledge about the unspecified risk factors. Both t and T are used to represent time. Using an extension of the notation in Table 1, values are needed for Fk(T,i) in order to simulate the polychotomous model as a function of time. However, the estimation method described in the preceding gives values for a smaller, dichotomous population. Thus, calling the latter cumulative probabilities q(T,i), it is desired to find the Fk(T,i). For each subpopulation, there will not only be the cumulative probability of an event, but also the cumulative probability of no event, rO,JT,i). To limit the size of the formulae in the following it is helpful to introduce the ratio pk(T,i) defined as Pk(T,i) =
ak(T,i)
1.0 - q(T,i)
(1)
Since all members of the subpopulation either experience the event or do not, the probabilities of these two alternatives must add to one. Symbolically this is, 4T,0
+ T~&T,~)= 1.0
236
E. Ackerman /ht.
J. Biomed. Compur. 36 (1994)
Pk(Ti)
This allows redefining the preceding ratio as Fk(T,i) Pk(T,i) =
‘RkvJ)
Fk(T,i) Fo(Ti)
or
(2)
Fd T,i) = pk(T,i) * FdT,i)
Note that the term Fo(T,i) is the cumulative probability of having no events and could have been written S( T,i). The sum of the cumulative probabilities of the polychotomous model must add to 1.0. In symbolic form, using m for the number of possible outcomes other than no event, k = m
Fk(T,i) = 1.0
(3)
k=O
Replacing the values for all but the no event probability in Eq. 3, by the equivalent expressions in Eq. 2, leads to k = n,
Fo(T,i) + c
[Fo(T,i) * &(T,i)]
= 1.0
&=I
Rearranging terms and solving, gives 1.0 Fo( T,i) =
(4)
k = ,,I
1.0 +
(5)
k = ,,I
&(T,i)
&=I
It is then assumed that this ratio also would have applied to the individuals who actually experienced other events. Therefore one may write
c
=
1.0 + c
rOk(r'i)
Pk(T,i) =
229-238
c
&(T,i)
&=I
Now it is an easy step to substitute the expression in Eq. 4 into Eq. 2 to obtain,
The preceding approximation, Eq. 5, is used repetitively throughout this series of tutorials. Its more complete justification and development can be found in Ref. 10. This type of approximation seems intuitively pleasing for CHD models where the F values are small for all of the events studied; other approaches would be needed if the F values were larger. In words, the approximation allows the estimation of the normalized risk coefficients using standard packaged programs, since the resultant expressions for the dichotomous risks can be suitably combined to find the estimated population-based risks. 8. SUMMERSand CRISPERS Once the required preliminaries have been completed, it is necessary to have a computer program if one wishes to simulate a micropopulation. Initially, the National Micropopulation Simulation Resource developed a series of models for specific applications and further specialized models within these areas. It became clear that a more general approach was needed to reduce redundant programming efforts and to allow greater efforts to be focussed on the simulation studies themselves. An abstract generalized model called SUMMERS was created for that purpose [5,11]. It and the associated SUMMERSSimulation Shell have proven successful over the years and have been used in numerous studies. One of the specializations of SUMMERS deals with CHD models; it is called CRISPERS. It formed the basic tool for two PhD theses [12,13] and for several journal articles [2,14- 181. Several additional studies are in progress using CRISPERS. The simulation programs based on SUMMERS allow the user to alter the nature of the model simulated. In the case of CRISPERS, the user is permitted to input information to alter the number of states and their interconnectivity. The user can also modify the population characteristics and the
E. Ackerman /In!. J. Biomed. Comput. 36 (1994) 229-238
parametric values associated with event probabilities. It is possible for the user to request several different types of interventions and a variety of reports at the end of the simulation. These features are realized because so much of the simulation is table driven, the values within the tables being easily modified at the time of program execution without the need to recompile. At the present time, the program has to be modified and recompiled if the probability rules are changed. Several projects are underway to remove this and other limitations as well as to further improve the user interactions. The advantages of using CRISPERSas a simulation tool are primarily those of convenience for the user. The latter is allowed to focus on the questions of research interest per se rather than on the details of programming. The user can also emphasize design of simulation studies confident that the programs used are efficiently written and carefully tested. The design, testing and validation of simulation programs can be quite time consuming, especially for the scientist who does not wish to become a full-time programmer. The SUMMERS simulation shell and the CRISPERS programs made it possible to carry out the simulations described in several of the references and also enabled the added studies now in progress. The CRISPERSprograms are referenced in the next five tutorials. However, they move away from the more general introductory remarks found here and into more substantive discussions. In particular, the next tutorial is concerned with the representations of the probability of the occurrence of events in models of interest in chronic disease epidemiology. The CRISPERS system is currently in use by several investigators at the University of Minnesota. Copies of the programs and sample runs are readily available from the National Micropopulation Simulation Resource.
author has discussed this paper with various colleagues, including Dr. La&l C. Gatewood, Ya-Jung Tsai, and Dr. David Gilbertson. The efforts of Jan Marie Lundgren in the preparation of the final manuscript are acknowledged. References I
2
3
4
6 7 8 9
10
II
12
13
14
Acknowledgements This work was supported in part by NIH Grant P41-RR01632. This paper includes information and descriptions of simulations which have formed parts of several papers and doctoral theses. The
237
15
Gail MH: A bibliography and comments on the use of statistical models in epidemiology in the 1980s. SIU/ Med. 10 (1991) 1819-1885. Zhuo Z. Ackerman E, Gatewood LC, Kottke T, Wu S and Park HA: Polychotomous multivariate models for coronary heart disease simulations. 1. Tests of a logistic model, In1 J Biomed Compur. 27 (1991) 133-148. Marmot M and Elliott P (Eds): Coronary Heart Disease Epidemiology: From Aetiology IO Public Health. Oxford University Press, Oxford, 1991. Ackerman E, Elveback LR and Fox JP: Simulation oflnfectious Disease Epidemics, C.C. Thomas. Springfield IL, 1984. Ackerman E, Zhuo Z, Altmann M. Kilis D, Yang JJ. Seaholm S and Gatewood L: Simulation of stochastic micropopulation models. 1. The SUMMERS simulation shell, Comput Biol Med. 23 (1993) 177-198. SAS: SAS UserS Guide; Staristics, 5th edn.. SAS Institute Inc. Cary NC, 1988. Dixon WJ et al. (Eds): BMDP Sratistical Software Manual, University of California Press, Los Angeles, 1988. IMSL: FORTRAN Subroutines for Mathematical Applicarions, IMSL Inc, Houston TX, 1989. Gross AJ and Clark VA: Survival Distributions: Reliability Applications in the Biomedical Sciences, John Wiley & Sons, New York, 1975. Begg CB and Gray R: Calculation of polychotomous logistic regression parameters using individualized regressions, Biomerrika, 71(l) (1984) 11-18. Seaholm SK: Adaptable Software for Simulation of Interacting Population Models, PhD Thesis, University of Minnesota, Minneapolis, 1987. Park HA: Simulation of a Population-Based Model of Coronary Heart Disease Morbidity and Mortality. PhD thesis, University of Minnesota, Minneapolis, 1987. Zhuo Z: Expert Support Sysiem to Select, Model and Interpret Cardiac Risk Funcfions, PhD thesis, University of Minnesota, Minneapolis, 1991. Zhuo Z, Ackerman E, Gatewood LC, Kottke T, Wu S and Park HA: Polychotomous multivariate models for coronary heart disease simulations. II. Comparisons of risk functions, Inf J Biomed Comput, 28 ( I99 I ) I8 I-204. Zhuo Z, Ackerman E. Gatewood LC and Kottke T: Polychotomous multivariate models for coronary heart disease simulations. III. Model sensitivities and risk factor interventions, Int J Biomed Compuc. 28 (1991) 205-220.
238
E. Ackerman /Ini. J. Biomed. Comput. 36 (1994) 229-238
16
Kottke TE, Gatewood L, Wu SC and Park HA: Preventing heart disease: Is treating the high risk sufftcient? J
17
Kottke TE, Gatewood L and Park HA: Using serum cholesterol to identify high risk and stimulate behavior change: Will it work? Ann Med, 21 (1989) 181-187.
Chin Epidemiol, 41 (1988) 1083-1093.
18 Zhuo Z, Ackerman E and Gatewood LC: An expert systern for simulation of coronary heart disease risk factor interventions. In Fifteenth Annual Symposium on Computer Applications in Medical Care (SCAMC) , American Medical lnformatics Association, Washington DC, 1991.
(Received 31 August 1993;accepted 8 September 1993)
Abstract
This is the first in a series of tutorials concerning the simulation of micropopulation models for use in epidemiological research. The series emphasizes the techniques that are in use at the National Micropopulation Simulation Resource at the University of Minnesota. This initial series of tutorials uses chronic disease models, specifically coronary heart disease models, to illustrate the principles and methodologies employed. The series is intended for persons having some familiarity with electronic, digital computation and desiring to know more about a specific class of simulations, namely those that apply Monte Carlo techniques in the simulation of micropopulation models as a part of biomedical research. It is the author’s hope that these tutorials will be particularly useful for students in medical informatics; they could also be of interest ,to professionals in both the computational and the biomedical sciences. Keywords: Simulation;
Micropopulation;
Risks; Stochastic model; Chronic disease model; Monte Carlo
1. Simulation Simulation is an extension of the concept of modeling complicated systems by simpler ones; this has been an important technique in the natural sciences for many years. For example, physicists and engineers work with mathematical representations of frictionless systems to learn about the behavior of real objects. Biologists work with animal models to represent human organ systems and * Corresponding author. E-mail: [email protected]. umn.edu.
with mechanical models to represent circulatory dynamics. The use of models in no way obviates the need for careful experimental observations; in fact, quite the opposite is true. Nonetheless, the models, by their very simplifications, make it easier for humans to understand natural phenomena. For many years, there was a widespread suspicion of mathematical models within biological communities. This distrust was frequently based on the belief that the simplifications necessary to arrive at a mathematically analysable model
0020-7101/94/~07.00 0 1994 Elsevier Science Ireland Ltd. All rights reserved SSDI 0020-7101(93)00939-F
230
E. Ackerman / Int. J. Biomed. Comput. 36 (1994) 229-238
overlooked too many key features of the real system. All too often that was a correct interpretation. However, since the advent of computer-based analyses and simulations, over-simplifications are no longer a necessary part of modeling. Computer-based simulation has gained respectability within the natural sciences and the social sciences. The applications of simulation techniques are so widespread, and the models so diverse, that it would not be appropriate to catalogue them in this tutorial. Neither would it be very instructive, since one could read a list of applications with no understanding of the simulations involved. Rather, greater detail is presented here about a limited group of simulation models. In some discussions of simulation, it is convenient to distinguish between deterministic models and stochastic ones. In the former, there is no uncertainty concerning the outcomes of the simulations, whereas in the latter the outcomes are distributions of events and of uncertainties about those distributions. In actual practice, the ideal of a purely deterministic or a purely stochastic model is never met; all actual models have some elements of both. This set of tutorials is concerned with models that are primarily stochastic in nature. The models presented are called micropopulation models in that each individual in the population is followed through a sequence of time periods. The models discussed are ones that are represented by difference equations rather than differential ones; the outcomes of these models are determined using Monte Carlo techniques. These terms are defined and examined more fully in this first tutorial. Although the restrictions in the preceding paragraph exclude most uses of simulation in biomedical research, the possible areas of the models described remains very broad. For pedagogic purposes, it seems suitable to limit such applications still further. Accordingly, this set of six tutorials describes the use of simulation in epidemiological studies of coronary heart disease. This first tutorial is meant primarily as an introduction to this type of simulation study and the last of the six tries to place it in a broader context. The four tutorials in between consider various aspects of these simulations. A frequent use of modeling in epidemiology is
concerned with macropopulation models that are primarily statistical in nature. These will not be discussed further here. The interested reader can find a comprehensive bibliography with comments in Ref. 1. The studies described in the current set of tutorials are on-going, and it is not possible to predict their future directions. Although emphasis is placed on completed studies, none the less studies in progress are mentioned. The results of the continuing studies should appear eventually in journal article format. Additional tutorials, however, may focus on other applications of the techniques employed at the National Micropopulation Simulation Resource. 2. Chronic disease -
coronary heart disease
Epidemiology originally meant a study of epidemic occurrences of disease. Epidemic, in turn, meant that it was outside of, or in excess of, the number of cases of that disease expected in the given population. Epidemics might be due to infectious agents or to transitory changes in the environment. Over the years the term epidemiology has been changed and extended to mean the study of the occurrence and etiology of disease in a population. The present series of tutorials is concerned with simulations of models of interest in epidemiology. In particular, this series has been restricted to models of chronic diseases, that is to ones that continue over an extended period of time. The diseases considered are non-infectious ones. For simplicity, only non-infectious diseases associated with the heart and the coronary circulation are considered. These are described by the label coronary heart disease and will be referred to by the letters CHD; the abbreviation is used in this fashion in numerous references, one of which is Ref. 2. A review of coronary heart disease epidemiology can be found in Ref. 3. For these simulation studies, it is convenient to consider a population as composed in part of persons who have not shown symptoms of CHD. Others in the population may have heartassociated pains called angina with or without any other evidence of CHD. The population will also be considered to include individuals who have
E. Ackerman /ht.
J. Biomed. Comput. 36 (1994)
experienced one or more non-fatal myocardial infarcts. The latter are occurrences in a region of the heart muscle myocardium in which the local circulation fails resulting in a permanent loss of function in that region of the myocardium. Where convenient, the abbreviation MI will be used to designate myocardial infarct. Other types of CHD are well known, but will not be considered in this series of tutorials. Rather, it will be assumed that all individuals are free from CHD except for possible angina and previous MIS. Some MIS are called silent because the individual does not report any symptoms; silent MIS are known only from autopsy data. In contrast to these, all of the MIS considered in the models presented are so-called acute MIS (AMIs). The latter may be fatal or non-fatal and may or may not involve hospitalization. Amongst the other events that may happen to a simulated individual are ones that involve death. Death may occur due to causes other than CHD; this type of event is referred to by the letters OCD. Death due to CHD (CHDD) may be preceded by, or the consequence of, an AM1 or simply occur suddenly with other etiology. The latter is referred to as sudden CHD death. The probability of experiencing CHD is highly age-dependent. Age may be specified as a risk, or the population simulated may be restricted to individuals within a given age range. In addition, at least in the past, the probability of CHD has been found to be quite different for men than for women. Thus, it may be necessary to carry out the simulations for sub-populations of only one sex. Heterosexual populations can be resolved into two separate ones for simulation purposes, or the simulation model may use the individual’s gender as a pointer to a sub-model. 3. Goals for coronary disease models The selection of a particular type of simulation model and of a set of simulation techniques is highly dependent upon the goals of the investigator. One goal that is emphasized in this set of tutorials is predicting the changes that might be anticipated as a result of population-based intervention programs. Another goal might be the evaluation of various risk factors, both in their
229-238
231
quantitative importance and also in the relative strength of their effects. A quite different approach is to seek a model that allows one to incorporate knowledge of the physiological processes underlying the pathologies modeled. By and large, in the area of chronic disease epidemiology, such models have not been widely used. By contrast, in infectious disease epidemiology, physiologically-based models are more prevalent [4]. However, CHD models that incorporate even a small approach to representing physiological processes are more intellectually satisfying than ones that make no attempt to do so. Another goal, albeit one that is often not mentioned explicitly, is to have a model that allows an easy estimate of the prevalence or the incidence of CHD-related events. The ease may depend upon the availability and user-friendliness of packaged statistical programs. Closely related to that goal is the one of simplicity. If two models are equally successful in other regards, the simpler one appears preferable. (The criterion of simplicity is sometimes called Occam’s razor.) Other goals are also important. One concerns the design of new surveys and new clinical trials. The use of a predictive model can help select the items to be determined, the population size and the duration of the surveillance. Similarly, model simulations can suggest the most critical parameters to monitor. In turn, the new studies can help to test the selected model. A frequent need is the communication of the results of a study to others. The use of models and of micropopulation simulations can make this communication easier. Thus the model and the simulation become a conceptual framework rather than a separate investigational tool. Goals will be referenced repeatedly in the discussions of specific models. By and large, a given model may be optimal for one of the goals cited in the preceding. Few, if any, models are well suited to all of the aforementioned goals. 4. Discretized stochastic micropopulation
models
In CHD studies, emphasis is placed on the occurrence of one or more of several possible events. Events, such as the occurrence of an MI, are dis-
E. Ackerman
232
/ ht.
J. Biomed.
Crete ones. Individuals experiencing such an event are said to move to a new state. Population models frequently involve simulations in which the fraction of the total population in a given state is considered. A model of that nature is said to be a macropopulation one. For such models, no attempt is made to keep track of an individual population member or an individual’s past history. On the other hand, when sufficient heterogeneity is present in the population, it is easier to follow each simulated individual. This type of model is called a micropopulation one. Only micropopulation models are considered in detail in this series of tutorials. The states of a model of that nature are illustrated in Fig. 1, which shows a flow diagram for a stochastic, micropopulation model that could be simulated using suitable techniques. The model is simulated in general for a limited period of time, for example 7 years. The length of the simulation is referred to as an epoch. During a simulation, all individuals who can experience some event, for example an acute MI, are followed either until they experience an event or until the end of the epoch, which ever occurs first. By and large, there is a low probability of a given type of event, so the largest group of simulated individuals
I
Death
from
Coronary Heart
Disease
No Previous History of Myocardial Infarct
Fig. I. Flow diagram of a model of coronary heart disease. This is an extremely simplified model of three possible events that may occur to a healthy individual. This and other simplified forms are used in several of the following tutorials.
Comput.
36 (1994)
229-238
will start with no experience of CHD. and continue that way throughout the epoch. In micropopulation models, it is possible to represent the probabilities of events as continuous functions of time which may be described by differential equations. However, in the presence of population heterogeneity, it is easier to deal only with discrete events and discrete units of time. That is true for all of the models to be treated in this and following tutorials. If time is discretized, then there is no meaning to trying to represent the events as occurring within a part of the basic time interval. Even though the program representing a micropopulation model may consider the individuals serially (that is, one after the other), it is not possible for events occurring during that interval to influence changes occurring during that interval. Thus the models are well described by difference equations; values of interest are determined by values at the previous time interval. The difference equations may be represented explicitly or may be described algorithmically as part of the simulation program. 5. Monte Carlo and other simulation techniques Models of the nature discussed here are readily simulated by Monte Carlo techniques in which pseudo-random numbers are used to decide if a given individual experiences a possible event in a discretized time period. Fig. 2 illustrates the fashion in which a uniform random variate can be used to select amongst several possible discretized outcomes. Effectively, this method transforms a uniformly distributed random variate to one that is distributed in the fashion appropriate to the simulation being considered. (It is also possible to transform a uniformly distributed random variate to any arbitrary continuous distribution using a technique similar to that illustrated in Fig. 2. Selection of a level or value based on a continuous probability function is used in some micropopulation simulations.) For illustrative purposes, assume that there are two possibilities from which one is to be selected. It is also necessary to assume that it is possible to compute the probability of each possibility. (The sum of the two probabilities must be 1.0.) Then a
E. Ackerman/Int.
J. Biomed. Comput. 36 (1994)
Outcome
I L-u t Event A
No
Change
0.0
0.2
Oi
Event B
0.6
Event C
08
1.0
Value of Random Variate
Fig. 2. Selecting outcomes using Monte Carlo techniques. The figure illustrates the process of choosing of one of four discrete options. First, a value is selected for a pseudo-random variate distributed uniformly in the range [O,l). If the value of the pseudo-random variate were 0. IS, i.e. below 0.4, the option No Change would have been selected, whilst if the value were above 0.8 the option Event C would have been chosen. Effectively this has transformed the uniformly distributed variate to a quadrunomial variate with probabilities 4 = 0.4, Pa = 0. I, P, = 0.3, and PC = 0.2.
random number uniformly distributed between 0.0 and 1.0 is selected. If the random number is smaller than the probability of the first possible outcome, that outcome is selected. Otherwise the second outcome is chosen. For more details please see Ref. 5. Where convenient, a wide variety of so-called packaged programs can be used to assist in the necessary statistical computations. Sets of packages that are referenced in successive tutorials include SAS [6], BMDP [7], and IMSL [S]. It is assumed that these are readily available to any group simulating stochastic micropopulation models representing the incidence and/or prevalence of CHD events in a given population or subset of the population. It is also assumed that any individual reading about such populations has at least an awareness of these specific packaged programs. The resulting model may be thought of as having at least three parts. One is comprised of the states and their interconnectivity as illustrated in Fig. I. This first part is sometimes referred to as a compartmental model since the states may be regarded as compartments of the population. However, the flow of individuals from one compartment to another is quite different than that of the usual meaning of compartmental analysis.
229-238
233
The second part of the CHD model is an algorithmic description of the procedures by which the probabilities implicit in Fig. 1 are determined. In some usages an algorithmic part is referred to as the statistical model, although this may be a quite misleading designator. A variety of algorithms are discussed in the next tutorial. These algorithms are described using one or another of the terms defined in Table 1; the notation used is modified from that found in Ref. 9. The third part of the CHD model of coronary heart disease presented is concerned with specific risk factors and their interactions. The specific risk factors included are chosen based on biomedical and epidemiological knowledge. The topics concerned with risk factors, their descriptions and their use in simulations are commented upon below in this tutorial and are developed more fully in following tutorials. 6. Field observations -
polychotomicity
The coronary disease models, since they are population-based and not driven by either known physiology or anatomy, must be verified by comparison with data on one or more real populations. The data are usually acquired by some type of survey or clinical trial. It is customary to refer to such observations as field data. Because of the fashion in which it is gathered, the data have far greater uncertainties than those associated with a laboratory experiment or with a controlled environment. The events of interest in coronary disease models occur relatively rarely. Thus the observations must be extended over a considerable period of time. Partly as a result of this, an individual is at risk for more than one possible type of event, several of which may exclude the others. This makes it harder to estimate the probabilities associated with each kind of event. Situations where only two outcomes are possible are called dichotomous, whilst those with multiple types of outcomes are designated polychotomous. One difficulty with all polychotomous models is estimating the parameters of the model. This is exacerbated by the mutual exclusivity of the events considered. One method of avoiding this limitation is to consider several subsets of the popula-
E. Ackerman I ht. J. Biomed. Comput. 36 (1994) 229-238
234
Table I Descriptors for dichotomous probabilities Name
Symbol
Event Density
Equivalent names
Definitions
Other relationships
event distribution pdfa death density
Pdr =f(t)dr Pdr is event probability during (1,~+ dt)
f(f) = dF(t)ldr f(‘,‘= -dS(r)/dt f(r)df s
s
= I
0
T
Cumulative Event Distribution
Survival Function
F(T)
S(T)
cdf” cumulative probability
survivorship death function
F(T) =
f(+ 0
-f(r)dr
S(T) =
s
T
F(T) = I - S(T) F(O)=0 F(a)
*
1
S(Q=l-F(T) S(0)= I S(m) * 0
Hazard Rateb
X(r)
force of mortality failure rate instantaneous death rate
f(l) = W)W) X(r) = -(dS(r)ldr]/S(r)
s
k(t) = -d[log(S(r))]ldr MO) =f(O)
T
Cumulative Hazard
A(T)
A(T) =
X(T)dr
0
A(r) = -logIS( T)l S(r) = e-*(n A(0) = 0 A(a) * m
“Abbreviations: pdf, probability density function; cdf, cumulative distribution function. bFor some purposes it may be desirable to use a slightly different definition for the hazard rate, namely A(T + 7,T) =/(T + 7)/S(T) which is the probability density of dying at time T + 7 conditional on survival to time T. This form is not used in these tutorials.
tion separately. Each smaller population contains the individuals experiencing no events or those individuals who experience a particular kind of event. To simulate the CHD model including the possibilities of different types of interventions, it is necessary to know not the probability of an event within the smaller population, but the probabilities within the entire population. Attempts to estimate these often fail for purely numeric reasons. Accordingly, an approximation is needed. The method described in the following section is valid only for events with a quite low probability. It permits estimations of the necessary parameters for polychotomous models, based on the more readily computed estimates from a group of dichotomous models. There are a number of different statistical terms that are used to describe the probability of an event. These are presented here for an individual in a dichotomous situation where the possibilities
are either an event that removes the individual from the population or no such event. These functions and distributions include probability density functions and cumulative distribution functions. The nature of these terms is such that if one is given, all of the others can be found, at least numerically. Table 1 presents a summary of these different functions; by and large, for the simulations of interest here, the cumulative distribution function (F) and the survivorship function (s) are the most important terms. A variety of expressions and rules can be used to describe the probability of the occurrence of an event in a given time interval, jI They are not commented upon further in this first tutorial, but are used throughout the series. 7. Multivariate risks factors
normalization of risk
It is sometimes pleasing to humans to assign a single cause to events that it is desirable to avoid.
E. Ackerman
/ht.
J. Biomed.
The cause is known as a risk factor. For coronary heart disease, it is necessary to include several different risk factors to simulate the events in the population. By and large, the risk factors have different units and are dimensionally incompatible. In order to combine them in mathematical expressions, it would be possible to multiply each risk factor by coefficients having the reciprocal dimensions. However, it then would be impossible to compare the relative importance of the risk factors since each coefficient would have different dimensions. In order to avoid that problem, one may use normalized dimensionless variates. A popular method of achieving this is to divide each risk factor by the population standard deviation. In order to further normalize the risk factors, it is also customary to subtract the population mean before dividing by the standard deviation. Thus a dimensionless, so-called normalized, risk factor represents the number of standard deviations that the actual risk factor is away from the population mean. This conversion may be expressed symbolically, letting yii represent the value of thejth risk factor for the ith individual. Then the corresponding dimensionless risk factor, x0 would be xv =
Yij - Yj OY,
where yj is the population mean for the jth risk factor and au. is its standard deviation. Although the unnormalized risk factors such as diastolic blood pressure and number of cigarettes smoked per day are always positive, the normalized values can be negative if the individual’s values are below the population mean. Dividing by some number such as the standard deviation is necessary to obtain dimensionless risk factors. However, the reference values need not be the population means; indeed, in some cases it makes more sense physiologically or epidemiologically if they are defined otherwise. The method of combining the risk factors to obtain an overall risk for a particular type of event such as an AMI is not dictated by any physiological or anatomic rule. Nonetheless, in the absence of any other indicator, it is customary for most models to use the exponential function of a linear combination of risk factors. This is a quite
Comput.
36 (1994)
229-238
235
arbitrary decision that is followed in most population studies. Provided the states and their interconnectivity, the algorithmic definition of the probabilities, and the form of the combination of the risk factors have been selected, it is usually possible to estimate the risk coefficients from the field data. As noted, this is most readily done by considering subpopulations with dichotomous outcomes. In the next few paragraphs, a method is described to estimate the risk coefficients for the entire population from the ones calculated for the subpopulations. In this tutorial, the normalized risk factors are represented by the symbol xii where the subscript i indicates an ordinal used to enumerate the individuals in the population and the subscript j, an ordinal used to enumerate the risk factors. Where appropriate the subscript k is used to indicate the particular outcome; the value k = 0 indicates no event. Since it is never possible to include all possible risk factors, the value j = 0 is used to indicate the combined influence of the unspecified risk factors. The value of xi0 = 1 is chosen arbitrarily for convenience in notation; it merely indicates a lack of knowledge about the unspecified risk factors. Both t and T are used to represent time. Using an extension of the notation in Table 1, values are needed for Fk(T,i) in order to simulate the polychotomous model as a function of time. However, the estimation method described in the preceding gives values for a smaller, dichotomous population. Thus, calling the latter cumulative probabilities q(T,i), it is desired to find the Fk(T,i). For each subpopulation, there will not only be the cumulative probability of an event, but also the cumulative probability of no event, rO,JT,i). To limit the size of the formulae in the following it is helpful to introduce the ratio pk(T,i) defined as Pk(T,i) =
ak(T,i)
1.0 - q(T,i)
(1)
Since all members of the subpopulation either experience the event or do not, the probabilities of these two alternatives must add to one. Symbolically this is, 4T,0
+ T~&T,~)= 1.0
236
E. Ackerman /ht.
J. Biomed. Compur. 36 (1994)
Pk(Ti)
This allows redefining the preceding ratio as Fk(T,i) Pk(T,i) =
‘RkvJ)
Fk(T,i) Fo(Ti)
or
(2)
Fd T,i) = pk(T,i) * FdT,i)
Note that the term Fo(T,i) is the cumulative probability of having no events and could have been written S( T,i). The sum of the cumulative probabilities of the polychotomous model must add to 1.0. In symbolic form, using m for the number of possible outcomes other than no event, k = m
Fk(T,i) = 1.0
(3)
k=O
Replacing the values for all but the no event probability in Eq. 3, by the equivalent expressions in Eq. 2, leads to k = n,
Fo(T,i) + c
[Fo(T,i) * &(T,i)]
= 1.0
&=I
Rearranging terms and solving, gives 1.0 Fo( T,i) =
(4)
k = ,,I
1.0 +
(5)
k = ,,I
&(T,i)
&=I
It is then assumed that this ratio also would have applied to the individuals who actually experienced other events. Therefore one may write
c
=
1.0 + c
rOk(r'i)
Pk(T,i) =
229-238
c
&(T,i)
&=I
Now it is an easy step to substitute the expression in Eq. 4 into Eq. 2 to obtain,
The preceding approximation, Eq. 5, is used repetitively throughout this series of tutorials. Its more complete justification and development can be found in Ref. 10. This type of approximation seems intuitively pleasing for CHD models where the F values are small for all of the events studied; other approaches would be needed if the F values were larger. In words, the approximation allows the estimation of the normalized risk coefficients using standard packaged programs, since the resultant expressions for the dichotomous risks can be suitably combined to find the estimated population-based risks. 8. SUMMERSand CRISPERS Once the required preliminaries have been completed, it is necessary to have a computer program if one wishes to simulate a micropopulation. Initially, the National Micropopulation Simulation Resource developed a series of models for specific applications and further specialized models within these areas. It became clear that a more general approach was needed to reduce redundant programming efforts and to allow greater efforts to be focussed on the simulation studies themselves. An abstract generalized model called SUMMERS was created for that purpose [5,11]. It and the associated SUMMERSSimulation Shell have proven successful over the years and have been used in numerous studies. One of the specializations of SUMMERS deals with CHD models; it is called CRISPERS. It formed the basic tool for two PhD theses [12,13] and for several journal articles [2,14- 181. Several additional studies are in progress using CRISPERS. The simulation programs based on SUMMERS allow the user to alter the nature of the model simulated. In the case of CRISPERS, the user is permitted to input information to alter the number of states and their interconnectivity. The user can also modify the population characteristics and the
E. Ackerman /In!. J. Biomed. Comput. 36 (1994) 229-238
parametric values associated with event probabilities. It is possible for the user to request several different types of interventions and a variety of reports at the end of the simulation. These features are realized because so much of the simulation is table driven, the values within the tables being easily modified at the time of program execution without the need to recompile. At the present time, the program has to be modified and recompiled if the probability rules are changed. Several projects are underway to remove this and other limitations as well as to further improve the user interactions. The advantages of using CRISPERSas a simulation tool are primarily those of convenience for the user. The latter is allowed to focus on the questions of research interest per se rather than on the details of programming. The user can also emphasize design of simulation studies confident that the programs used are efficiently written and carefully tested. The design, testing and validation of simulation programs can be quite time consuming, especially for the scientist who does not wish to become a full-time programmer. The SUMMERS simulation shell and the CRISPERS programs made it possible to carry out the simulations described in several of the references and also enabled the added studies now in progress. The CRISPERSprograms are referenced in the next five tutorials. However, they move away from the more general introductory remarks found here and into more substantive discussions. In particular, the next tutorial is concerned with the representations of the probability of the occurrence of events in models of interest in chronic disease epidemiology. The CRISPERS system is currently in use by several investigators at the University of Minnesota. Copies of the programs and sample runs are readily available from the National Micropopulation Simulation Resource.
author has discussed this paper with various colleagues, including Dr. La&l C. Gatewood, Ya-Jung Tsai, and Dr. David Gilbertson. The efforts of Jan Marie Lundgren in the preparation of the final manuscript are acknowledged. References I
2
3
4
6 7 8 9
10
II
12
13
14
Acknowledgements This work was supported in part by NIH Grant P41-RR01632. This paper includes information and descriptions of simulations which have formed parts of several papers and doctoral theses. The
237
15
Gail MH: A bibliography and comments on the use of statistical models in epidemiology in the 1980s. SIU/ Med. 10 (1991) 1819-1885. Zhuo Z. Ackerman E, Gatewood LC, Kottke T, Wu S and Park HA: Polychotomous multivariate models for coronary heart disease simulations. 1. Tests of a logistic model, In1 J Biomed Compur. 27 (1991) 133-148. Marmot M and Elliott P (Eds): Coronary Heart Disease Epidemiology: From Aetiology IO Public Health. Oxford University Press, Oxford, 1991. Ackerman E, Elveback LR and Fox JP: Simulation oflnfectious Disease Epidemics, C.C. Thomas. Springfield IL, 1984. Ackerman E, Zhuo Z, Altmann M. Kilis D, Yang JJ. Seaholm S and Gatewood L: Simulation of stochastic micropopulation models. 1. The SUMMERS simulation shell, Comput Biol Med. 23 (1993) 177-198. SAS: SAS UserS Guide; Staristics, 5th edn.. SAS Institute Inc. Cary NC, 1988. Dixon WJ et al. (Eds): BMDP Sratistical Software Manual, University of California Press, Los Angeles, 1988. IMSL: FORTRAN Subroutines for Mathematical Applicarions, IMSL Inc, Houston TX, 1989. Gross AJ and Clark VA: Survival Distributions: Reliability Applications in the Biomedical Sciences, John Wiley & Sons, New York, 1975. Begg CB and Gray R: Calculation of polychotomous logistic regression parameters using individualized regressions, Biomerrika, 71(l) (1984) 11-18. Seaholm SK: Adaptable Software for Simulation of Interacting Population Models, PhD Thesis, University of Minnesota, Minneapolis, 1987. Park HA: Simulation of a Population-Based Model of Coronary Heart Disease Morbidity and Mortality. PhD thesis, University of Minnesota, Minneapolis, 1987. Zhuo Z: Expert Support Sysiem to Select, Model and Interpret Cardiac Risk Funcfions, PhD thesis, University of Minnesota, Minneapolis, 1991. Zhuo Z, Ackerman E, Gatewood LC, Kottke T, Wu S and Park HA: Polychotomous multivariate models for coronary heart disease simulations. II. Comparisons of risk functions, Inf J Biomed Comput, 28 ( I99 I ) I8 I-204. Zhuo Z, Ackerman E. Gatewood LC and Kottke T: Polychotomous multivariate models for coronary heart disease simulations. III. Model sensitivities and risk factor interventions, Int J Biomed Compuc. 28 (1991) 205-220.
238
E. Ackerman /Ini. J. Biomed. Comput. 36 (1994) 229-238
16
Kottke TE, Gatewood L, Wu SC and Park HA: Preventing heart disease: Is treating the high risk sufftcient? J
17
Kottke TE, Gatewood L and Park HA: Using serum cholesterol to identify high risk and stimulate behavior change: Will it work? Ann Med, 21 (1989) 181-187.
Chin Epidemiol, 41 (1988) 1083-1093.
18 Zhuo Z, Ackerman E and Gatewood LC: An expert systern for simulation of coronary heart disease risk factor interventions. In Fifteenth Annual Symposium on Computer Applications in Medical Care (SCAMC) , American Medical lnformatics Association, Washington DC, 1991.