We begin our exploration of continuous, deterministic, compartmental models in epidemiology with the simplest model possible. We examine a model of a fixed population with no births or deaths or migrations that experiences a constant rate of disease. This enables us to expound upon the relationships between risks and rates while clarifying some of the most basic abstractions and methods of the modeling system we will employ. It also provides an introduction to the use of Stella™. The model is simple enough so that there is a "closed form" solution to the differential equations behind it, allowing us to elaborate upon the relationships between closed form models of time effects and dynamic simulations.
Just as one can successfully execute the mechanics of analyzing data without fully understanding all the abstractions involved and their implications, the same is true for building models. The reason you will build models in this course, however, is to better comprehend of epidemiological concepts so that you have a sound basis for constructing epidemiological theory. We are pursuing a very hands on approach to model building in which we postpone much of the theoretical and philosophical discussion about our endeavor until we have some broader experience with models. Before jumping into the mechanics, however, we will ponder some of the most fundamental concepts underlying the compartmental models we will construct. These have to do with the time distribution of disease as we conceptualize them as risk or rates.
Risks and rates are theoretical abstractions (models) with which we try to perceive and make sense of disease patterns in populations. The nature of these abstractions and the simplifying assumptions we make when using them are fundamental to the practice of epidemiology.
Our model of risks employs employs the concept of dimensionless proportions whose values range from zero to one. Just as probabilities may be interpreted in relative frequency or personal probability terms, risks can have population frequency or personalist interpretations (See for example Statistical Inference by Michael Oakes, Epidemiology Resources Inc. 1990). In the personal conceptualization, risk represents the danger that something will happen to an individual. In the frequentist interpretation, disease risks are the fraction of times that disease occurs. If 20% of at risk individuals present at the start experience disease during a particular interval of time, then the risk of disease is 20%. The frequentist calculation of risk is often translated into a personal risk probability by assuming that an individual is representative of the group in which the frequentist risk was calculated. The essential abstraction in conceptualizing risks as they are estimated and used in epidemiology is the assignment of individuals to groups.
Although risks are mathematically dimensionless, they have an implicit time period over which they occur. Epidemiologists often fall into the bad habit of ignoring this implicit time period. In a practical sense, this time period provides a dimension to risks, even though the mathematical expression of the risk is dimensionless.
Our model of rates employs the concept of continuous processes that theoretically is going at some speed at any point in time. Rates, however, are never estimable at any point in time in the same way that we can look at the speedometer of our car to get an estimate of its speed. They are estimated by counting the number of events per person years of observation of individuals at risk of the event. They have the dimensions of events (case onsets) per person-years and may have any value from zero to infinity. When is a rate of 1,558,000,000,000,000 the same as a rate of 1? When the rate of one per million person centuries is expressed in terms of person seconds. (1.558 times 1015 person seconds equals one million person centuries). Even though the numerical value of a rate may change dramatically as we change the denominator units, the rate does not change. You see, rates are not dimensionless like risks are.
Whereas the calculation of risks requires the same abstractions as the practical uses of risks, the way rates are calculated does not require the same abstractions as the way we will use them in constructing models. To calculate rates, we merely count events in the same way we count events to calculate risks. The way we will use rates in the models we will construct, however, the events in the rates are not composed of countable entities. They are continuous rather than discrete.
We will often think about rates in the models we will build in terms of what is happening to individuals. But the models we will build will not be of individuals. They will be of abstract populations in which no individual distinctions are made. We do not have time in this course to cover models with discrete individual units. Our modeling technology will make the abstraction that our population is infinitely divisible. We will model populations as separate compartments. Population will flow between compartments. Formally we should only refer to the size of the "compartments" or "stocks" which we will model as sizes of populations. But we will think of individuals flowing between compartments rather than some abstract continuous entity called a population. Thus we will find it useful to treat the size of stocks or compartments as number of individuals. But our modeling approach assumes that the population is infinitely divisible. Therefore when we refer to the size of a compartment as number of individuals, we will generally have fractions of individuals. We will have to learn when this somewhat unrealistic abstraction will distort the phenomenon in which we have a scientific interest and when we can use its power fully. Sometimes population phenomena such as temporal patterns of person to person spread epidemics can be quite different for models with discrete individuals than for continuous models like those we will use in this course. For most chronic disease phenomenon, however, continuous and discrete model results will not differ so much.
The way epidemiologists most commonly use the concept of rate is to assume that it is constant over time. This assumption, together with the way epidemiologists estimate risks by counting events over discrete time intervals, makes the use of rates much more like the use of risk in epidemiological thinking. In using compartmental models to construct epidemiological theory, however, we use several abstractions that make rates behave quite differently from risks.
Most introductory text books in epidemiology do not go into the relationships between rates and risks. Some do not even make any distinction between them. They develop all of the pertinent epidemiologic measures in terms of risks and disregard rates, even though they may use the word rate when referring to some of their risk concepts. Thus many epidemiology students may only relate to disease development at a population level in terms of risk. The concept of an instantaneous rate in a population is elaborated only in more advanced texts such as Epidemiologic Research by Kleinbaum, Kupper, and Morgenstern, or Modern Epidemiology by Rothman.
Perhaps one reason for the early focus of epidemiology on risks rather than rates is that risk calculations are more widely applicable than rate calculations. Risk calculations are meaningful for point source and person to person spread epidemics while rate calculations are often not relevant in these situations. The risk of getting ill during an epidemic can be readily calculated and presumably will be unchanging once the risk has passed. In contrast, the rate of illness in an epidemic, be it a point source epidemic or a person to person transmitted one, usually rises from zero to a peak and then falls to zero again. There is no single rate that can characterize the disease process in an epidemic. Thus the presentation of risks is more comprehensible in these situations than the presentation of rates.
A major reason risk concepts dominate over rate concepts in epidemiology is that risk concepts underlie almost all of the theoretical advances that have been made in epidemiology. Risks are probabilities where each of the units experiencing an event in the numerator is also found in the denominator and only units with a conceptual chance of getting into the numerator are included in the denominator. Most of the statistics developed to deal with epidemiological data are based upon probability concepts and theory. Even more importantly, most of the fundamental concepts in epidemiology such as confounding, effect modification, information bias, sample selection bias, attributable risks, and population attributable risks, are founded upon risk calculations and probabilities.
One reason the concept of risk is more natural to epidemiologists than the concept of rate is that risk deals with discrete events occurring over discrete time intervals. The concept of rate, on the other hand, requires continuity abstractions to which it is not easy to relate and which cannot be manipulated productively without a background in calculus. The numerator and denominator in rates are assummed to be continuous entities. That is why we can translate a rate of one per person per year into a rate of 1/365 per person per day. The essential abstraction in the concept of rate is that instead of thinking about what happens over an interval, we think about the state of affairs at a point in time.
We can think about something happening over a short time interval dt. To make the abstraction to what is happening at a point in time with no duration, we have to make the abstraction that we keep making dt shorter and shorter until it finally reaches the limit where it has no duration. Conceiving of what is happening at a point in time with no duration requires a conceptualization of time as being continuous. That is to say a time interval can always be further divided, no matter how short it is. Nothing can actually happen at a point in time that has no duration. So what does it mean to think about the rate at which it is happening as dt goes to zero?
To answer this question requires another abstraction which is even more strange--namely the continuity of events. The continuity of events implies that events can always be further divided. This abstraction is implicit to the concept of instantaneous rate. Since the event of disease developing occurs in individuals at separated points in time, the concept of continuous events implies that there are an infinite number of individuals available in whom events might occur. Assuming there is an infinite number of individuals is the same as assumming that the compartments or populations in our models are infinite. Such assumptions often seem quite unnatural to epidemiologists whose major concerns deal with data collection from individiuals and its analysis.
The reasons for the dominance of risk concepts in epidemiology seem clear. Why then do we concentrate on the concept of rate in this course? Why do we seek to alter a tradition that has been so productive using concepts to which our peers do not easily relate?
We believe that the science of epidemiology can be made even more productive and exciting than it is if we seek to develop a basic understanding not just of disease risks, but of the processes which generate disease in populations. We believe that it should be the major task of epidemiology to generate scientifically sound theory about the processes that generate patterns of disease in populations. We emphasize processes and the calculation of parameters which are scientifically generalizable on the basis of sound theory as an essential supplement to the more commonly pursued task of estimating risks which are generalizable only on the basis of statistical theory, which means they are generalizable only to the statistical target population from which the data was derived.
Processes imply dynamics or change. Rates are an underlying measure of process, of continual change. They provide a basis for constructing and analyzing models of processes which avoid many of the complications that arise when models deal with discrete individuals. If we used a risk over one time interval of 5 per thousand to say what happens to a population of 100 discrete individuals during that time interval, we would predict that 0.5 individuals would experience the disease outcome. We would either have to set up some kind of rule as to whether there were in fact no disease events or one disease event in these cases or we would have to bring chance into our models and use probability density functions to describe the outcome. In either case, this makes the models very very difficult to deal with. The abstraction of a continuous rate simplifies things greatly.
Whereas the concept of risk may be more widely used in epidemiology than the concept of rate, the concept of rate is theoretically more fundamental to our models because in them rates generate risks and not vice versa. In the real world, it may be that risks experienced by individuals generate rates as we calculate them. But in our models, rates will generate risks. A rate acts at an instant of time. As time moves across a time period, the cumulative action of the rate generates a risk across that time period. Using the concept of rate we can build more relevant and complete causal theory than we can using the concept of risk. We do this mathematically by making rates the measures upon which the integral calculus operates across time. Because risk is cumulative, it lumps different causal effects acting within a time period. Rates, on the other hand, are instantaneously responding to causal influences and more directly reflect the effects of individual causal factors.
One may question the theoretical utility of continuity abstractions intrinsic to the concept of rate because they seem so far from reality. But as one sees all the wonderful mathematical relationships derived making this sort of abstraction in differential and integral calculus and differential equations, the nature of the abstraction becomes natural and powerful. In this course you will not see many of these intrinsically beautiful and powerful relationships mathematically. We do not present much mathematics. Instead the abstraction intrinsic to calculus and differential equations will acquire a natural utility and meaning for you as we see the behavior of system simulations where we reduce dt until our system simulation behaves just as a truly continuous system would.
We hope that by being able to generate systems of continuous causes
easily, you will be better able to relate to cause as a system
concept. We hope this will make you much better in generating
hypotheses and theory about how patterns of exposures in populations
and the disease processes to which they relate generate patterns
of disease in populations.
Let us now do some practical model implementation of a constant rate process so that we can see how a rate generates a risk. Our purpose in this model construction is merely concept clarification. You might think that you already know all there is to know about risks and rates. But this model construction and analysis will serve to deepen your understanding in some respects. It will also help clarify what we mean by a compartmental model and how we use continuous compartments and continuous events to model conceptually discrete events occurring in discrete individuals. Another didactic purpose is that it will help clarify how we use discrete time intervals to approximate what we conceptualize as continuous time. And it will teach us a little about how Stella® works.
After familiarizing yourself with the operation of Stella™ using the material provided with the software, construct the following model at the "formula-diagram" level.
Stella™ stocks in our model will constitute compartments consisting of population segments. For this model we use "resevoir" stocks. Thus our populations are evenly mixing and there will be no distinction between any units in the population. Any unit in the well population could be the next unit to flow out into the diseased population but it will make no difference which unit is the next one to flow out. We will inevitably think of the units in the compartments as individuals but we will treat these population units as continuous as discussed earlier. The flow in our model is unidirectional. The well population can move into the diseased compartment but in this model there is no recovery from disease so the diseased cannot move back into the well.
What are the variables and what are the parameters in this model? The variables are model elements that change over time according to the processes or flows that we establish in our model. The parameters are elements that we fix. In the compartmental models of populations which we construct, the stocks or population compartments will always be variables. The value of those variables will be the size of the compartment. Note that this concept of variable is quite a bit different from the concept of variable as a measurement made on individuals. It is a population rather than an individual concept.
The variables in this model are the Well and Diseased compartments or stocks. For simplification purposes we now cut these down to W and D. This model has a single parameter. It is the rate at which disease develops in the well population. We will assume this rate is constant over time and give it a numerical value of 0.005 per person per time unit. Note that if we are using observational data to fix our rate parameter, we must be very careful to set our time units in correspondence to the data. The time units must be the same for every part of the model. One cannot use rates per month in one part of the model and rates per day in another and still have the model make sense. To simplify we now designate the rate parameter with just an "r".
If the compartments are variables and the rate is a parameter, what is the flow between well and diseased? The flows in compartmental models are merely functions of the parameters and the compartments. The modeler must define the form of these functions. The art of compartmental modeling of populations has several elements. The first is choosing how to divide the population into different compartments. The second is arranging the flows between compartments. The third is determining what parameters and what variable values affect the flows. This third element consists of placing connectors from compartments and parameters to the flow regulators. A final crucial element is formulating the relationships between variables and parameters that determine the flows. In our simple example this is done wholly within the flow regulator. Later we will have more complex formulations where we will put part of the relationships into Stella™ converters.
For simplification purposes in our model we will designate the flow from well to diseased as ƒW_D and we will formulate it as the size of the population on which the rate acts times the rate value. We use the asterisk for the times symbol so our flow is r*W.
The model we have constructed is of a closed population. The entire population is divided into two compartments and there are no flows into or out of the population. When there are such flows, Stella uses clouds to symbolize the source of population which flows into the population as well as the destiny of the population which flows out.
Note that we used a Stella™ converter for our model parameter. Converters may have several other functions in our models. We just mentioned that we may want to use them to construct relationships that are part of our flow formulations. Very commonly we will use them to define derived variables. A derived variable of interest in our model might is the disease prevalence or the fraction of the population which is diseased. We put a converter on our model and use connectors from W and D to this converter to make such a derived variable. Let us call this prevD. In this converter we formulate the conversion as D/(W+D). We then have the following diagram for our model:
There are no question marks on the derived variable because we already entered the formula D/(W+D) so that converter has now been defined. Likewise we defined the flow as r*W. So we have defined the relationships of interest in our model. Now we must define the remaining elements.
A decision to be made now is what size are we going to designate for our population. Is our population to have ten billion individuals or only ten? That might seem at first like an important decision. But in fact it is not because the population we are modeling is infinitely divisible. Whether we put the value of 10,000,000,000 or 10 on the overall population, the assumptions we will be using make the true number of individuals infinite. In the case of designating a total of ten billion for the population we might have 5 individuals in the D compartment and the rest in the W. In the case of designating the population size as 10, however, we might have 0.0000000005 "individuals" in the D compartment and the rest in the W. The proportion is the same. The fact that we can have 0.0000000005 "individuals" in a compartment speaks to the arbitrariness of the population size unit. We can convert population sizes merely by a process of division or multiplication. In our simple model here, the conversion is very simple but in complicated models it may be trickier. One must always be very careful to insure that the units in all elements of the model correspond to each other and in complicated models, it might be easy to leave some conversion undone or improperly done. It is for that reason, and that reason only, that the choice of initial population size is important. Choose a population size that will be meaningful to the audience to whom you wish to communicate. Then be sure to use units that are consistent in different elements of the model. We choose an initial population size of 10,000 and at the start we put all the population in the well category. This would be like starting a cohort study where initially only the well population enters the study.
For the general conceptual purposes that motivate this model construction, it does not make any difference what value we assign to the rate parameter. We will give it a value of one tenth per person per time unit. We will designate the time units to years. Note that the time units we designate does not change any values in our model nor the shape of any graphs we generate. The only time the choice of time units makes a difference is when we are using real parameters in our models for the purpose of projecting the behavior of real systems. Most often in this course we will be using our models more to improve our understanding than to project real world events.
While the value we assign to the rate parameter may be arbitrary, it does influence the value of the time step needed to generate model dynamics over time. In order for the numerical integration of our model to give values that correspond to those of the true differential equation, the fraction of population that flows out of a compartment across any time step must be small. If our rates have very high values and generate very high flows, then for the fraction of population flowing across a time step to be small, the time step must be very small. Stella does not allow for the value of the time steps to be too small. But the actual time steps can be made smaller by changing the time units rather than the value of the time step. At the same time step of 0.1 time units, the actual time step in month units will be one twelfth the size of a 0.1 time unit step in year units. Just because we may have calculated a rate from real data in year units does not mean we have to use the value in year units. We can just divide the rate by 12 and call it month units.
We have now completed our model. Let us examine the difference equations that Stella™ has generated from the model.
D(t) = D(t - dt) + (ƒW_D) * dt
INIT D = 0
INFLOWS:
ƒW_D = r*W
W(t) = W(t - dt) + (- ƒW_D) * dt
INIT W = 10000
OUTFLOWS:
ƒW_D = r*W
prevD = D/(W+D)
r = 0.1
Stella™ first generates equations for each stock with inflows and outflows. The value of the stock after the current time step is the value before the time step plus the inflows minus the outflows. Stella™ always formulates these inflows and outflows separately rather than substituting them directly into the difference equation. After the difference equations, Stella™ presents the converter formulas. Remember that some converters represent model parameters and other converters may just be functions of stocks and other converters.
In the model we have created Stella generates two difference equations, one for each stock. To describe our model we in fact need only one difference equation since the value of the other stock is defined by the formula W+D=10,000. Instead of using one difference equation and one algebraic equation to define our model, however, let us write both difference equations in a more common form. We substitute the values of the flows into the difference equations and the values of the parameters into the flows. Our model is then fully defined by the two difference equations
.........................................Equation
set (1)
Note that in these equations, the flow out of W equals the flow into D. That is consistent with the diagram we made of the model.
To numerically solve a set of equations like (1), a set of initial values at t-dt need to be established. Solving in this case refers to determining the Stock values and derived variable values over time. Once it has starting values Stella ™ uses one of three different numerical integration methods to "solve" the equations. We commonly present the solutions in terms of graphs or tables.
To convert the difference equation (1) into a differential equation, subtract the stock value in the prior time step from both sides and divide by dt. This gives the following
The numerater on the left is the change in W over the interval dt. As we let dt tend to zero, at the limit of dt tending to zero we form the differential equation
................................Equation
set (2)
The difference between equation sets (1) and (2) are that in (1), time is conceptualized as occurring in discrete steps while in (2) it is conceptualized as being continuous. The left hand side of the equations in (2) are read as "the instantaneous rate of change in W (or D) per unit change in time".
Stella™ always writes out its equations in difference equation form because the equations do not have to be used as differential equations as we will use them in this class. For many purposes the difference equations can be used directly. But the literature on epidemiological processes almost always presents differential equations. You should be able to write out the equations of the models you construct in standard differential equation format for presentation in the literature By the end of this class you should also be adept at taking a set of differential equations from the literature and implementing them and exploring them on Stella.
For practical purposes we cannot make dt infinitely small and still numerically solve Stella equations over time. How then can we use Stella in a differential equation mode? We make dt small enough for all practical purposes. We determine that dt is small enough for practical purposes by checking to see that the solutions of the difference equations do not differ as dt is cut in half. That is best done by putting data into a table at some broad intervals, saving the values from that table, redoing the simulation with the dt cut in half, and checking to see at various time points that stock values do not differ. If they do not differ, we can treat our difference equation solutions as differential equation solutions.
If they do differ, there are two things we could do. We could go to a higher order numerical integration procedure. The order in Stella from low to higher is Euler, Runge-Kutte 2, and Runge-Kutte 4. Then we could check at this higher numerical integration method to see if differences persist when the dt is cut in half. If they do, then further reduction of the dt is needed.
Solving or running a model produces output which defnes the values of compartments and derived variables at different points in time. Model output may be observed in graphs, tables, or animations. You should become familiar with all three by consulting the Stella™ manual. For the most part we will use graphs and tables. Stella™ produces this output by taking steps of size dt starting with the initial conditions. The value it gets at one point it then uses as starting conditions to project another dt step. When the Euler method is used, the slope used to project the change in the stock values is just the value of the differential equation given the starting stock values. If the slope is changing as the stock values change (which it almost always is) and if the steps taken are too large, then the difference equation solution calculated by Stella™ will be quite different from the differential equation solution. If we want to use the model as a numerical approximation to differential equations, we must set dt at appropriate values. Consider the following graph of an Euler solution where the curve is meant to represent the true differential equation solution and the linear steps represent an Euler approximation to this. If the step is too large for the rate of change of the slope, the difference equation solution will be quite different from the differential equation solution. Using the Runge-Kutta methods takes into account the rate of change of the slope. The Runge-Kutta 4 does this better than the Runge-Kutta 2 but it takes more computer time. The Runge-Kutta methods do not use the slope at the start of the interval but rather use a slope which is a type of average slope across the interval. These methods don't generate as much deviation as seen in the following graph. We can get the differential equation solution to better approximate the differential equation solution by using either a smaller dt or a higher order numerical estimation method.
Across each small dt in difference equations, the rate of a process per dt is treated as if it were the risk of the process. Using the Euler method it is the instantaneous risk at the start of the interval that is used as a risk. As the dt is shortened, the risk across the dt becomes less and less. There is a general rule of thumb that as the risk across an interval reaches 0.04, there will not be much difference of that risk from the rate per dt. NOTE WELL It is the risk approaching the value of 0.04 and not the rate which is involved in this rule of thumb. As we have seen earlier, any rate can assume any numerical value just by changing the time scale to which it refers. Here we want the time scale to correspond to DT.
Set up your model to run over 30 years using a dt of 0.125 years from its intial conditions and then set up a graph with the minimum and maximum values of D and W from 0 to 10,000, the min and max for PrevD from 0 to 1, and the min and max for the flow from 0 to 1000. Before running the graph Before running the graph Before running the graph draw in the shape of the curve that you expect the model to generate. A graph is provided here to help you in this task. Please do not be a slouch and run the model before you draw out what you expect. Consider whether the generated curves will be monotonic, if they will be monotonic, will they be convex or concave, and how fast and how far they will rise or fall.
to be FTPd to Dr. Koopman:
Please describe any difference between what you predicted the model would produce and what it actually produced. More importantly, explain why the model produced the curves that it did and why your predictions were not completely accurate.
Run the model at a dt of 0.5 years, 0.25 years, 0.125 years and 0.0625 using the Euler method and the Runge-Kutta 2 and Runge-Kutta 4 methods. Record the same variable values in the graph in a table at 15 year intervals so you have only three lines of results. For each method, save your table results in a spread sheet or word processor so you can compare them. Note the time the program takes to solve (run) the model. Try comparing instantaneous or summed table presentations with either beginning balances or ending balances. Determine which is the maximum time interval or minimal numerical solution method that you can use with this model and still have it perform for all practical purposes as a differential equation. To increase your understanding of how the numerical solutions work, uncheck the non-negative value box in the stocks, set the method to Euler, set the dt to 1, and increase the r to 0.5, then solve. Next try values for r of 1.0 then 1.5. Then do the same using the Runge-Kutta methods.
What is the most efficient way to solve this model so that it behaves as a differential equation model?
Identify the values in the three line table results that differ for different dts and solution methods. Explain why the differences occur without going into any analysis of how the solution methods work.
When differential equations are treated symbolically and not just as limiting situations for difference equations, they have a powerful base of mathematical theory which has been developed since the time of Sir Isaac Newton. It is this theoretical base that has made differential equation modeling attractive to mathematicians. It will be the exceptional epidemiologist who will deal with symbolic solutions of differential equations. Most epidemiologists will deal just with numerical solutions like those effected step by step using the Stella program. But it is useful for the epidemiologist to be aware of what the mathematician might have to offer that could advance epidemiological theory.
Difference equations can be examined theoretically to give insight into the dynamics of a system. Equilibrium points can be found. The local stability and sometimes the global stability of those points can be examined. By local stability we mean if the system gets pushed off the equilibrium by only a little bit, will it return to the equilibrium point. By global stability we mean will the system return to that point no matter how far it is knocked off from it.
Another advantage of differential equations is that sometimes there is enough mathematical theory available to find "closed form" solutions for them. With a "closed form" solution, if one is given the initial compartment values and parameters for the system, one can calculate directly what the compartment values will be at any specified time in the future without having to go through the small dt step by step numerical integration to calculate the stock values at any time in the future.
Our simple system of a constant disease rate in a fixed population has a closed form solution that relates risks and rates. When we start out with everyone in the well state in our model, the number of diseased individuals at any point in time over the total number of initial individuals is the cumulative disease risk up to that point in time. If we start out with only one individual in the well category, then we don't have to worry about dividing by the original population in order to get a risk.
From the theoretical arguments of calculus one can then solve the differential equations with the cumulative risk of disease on the left as:
.......................Equation (3)
"e" in equation (3) is the base of the natural logarithms. Instead of using the arguments of calculus directly to demonstrate this solution, we will use related arguments that demonstrate what the base of the natural logarithms actually is. For those students who have a very solid background in calculus, they may get more meaning out of the calculus arguments. But for those of us who are not powerful mathematicians, the following demonstrations provides insight into equation (3).
We begin by considering a compound interest problem. Suppose you
live in a country where the annual interest rate is 100%. If the
interest is compounded once a year, at the end of a year you will
have twice as much money as you originally started with. Let us
set up the problem this way.
M(12 months) = M(time 0) + 100%(M(time 0) = M(time 0) + 1(M(time
0) = 2(M(time 0)).
Now let us consider how much money we will have after 12 months
if the annual rate of interest stays at 100% but it is now compounded
every six months. If the annual rate of interest is 100%, the
six month rate will be 50% or 0.5.
M(6 months) = M(time 0) + 0.5(M(time 0) = 1.5(M(time 0))
M(12 months) = M(6 months) + .05M(6 months) = 1.5(M(6 months))
= 1.5(1.5(M(time 0))) = 2.25(M(time 0))
Now let us consider how much money we will have after 12 months
if the annual rate of interest stays at 100% but it is now compounded
every month. Analogous to what we have done above:
M(12 months) = M(11 months) + (1+1/12)M(11 months) = (1+1/12)(1+1/12)M(10
months) = etc. ...
= (1+1/12)12M(time 0) = 2.613M(time 0)
We could progress to weekly compounding, daily compounding ever
12 hour compounding, every six hour compounding, and instantaneous
compounding as per the series below.
n = 52 --> 2.693
n = 365 --> 2.715
n = 730 --> 2.71642233
n = 1460 --> 2.71735149
n = infinity --> 2.71828183... = e, the base of the natural
logarithms.
Note that as you double the number of steps, the differences between
the total value of the compounded funds becomes less and less
at higher numbers of compounding intervals. At some point when
you double the time interval you get no noticeable change at all
in the total. That is when we say the result has approached its
limit. In this case we give the limit the symbol "e".
The number "e" is a mathematical convention with which
you should become familiar called the base of the natural logarithms.
It is the base from which natural growth processes proceed when
instantaneous rates of growth are proportional to original size.
"e" is the base for both positive or negative growth.
Note that calculating interest across a single time step is like
calculating disease occurrence given a risk. Calculating interest
across any time interval in n steps with n = infinity is like
calculating disease occurrence in that time interval given the
rate which is effective in that interval.
In formal mathematical terms we can represent what we have done in the last example where we examined a series of increasing number of times when 100% interest was compounded over a defined interval as follows:
...............................................Equation
(4)
n is the number of times we compounded interest during the single
time interval across which we are calculating the accumulated
value of our savings.
A plot of n as a function of n follows:
The larger n becomes, the smaller is each time interval at which
we calculate accumulated interest and compound it. As the time
interval becomes very small the process is almost instantaneous.
In general in this course, we will think of rates as almost instantaneous
processes. The limit as "n" tends to infinity would
take us to a truly instantaneous process. That is the theoretical
world of calculus. All the simulations we will develop in this
course will be at very small intervals. We will in general try
to check to make sure that our intervals are so small that making
them any smaller will not change the result. This will make these
results equivalent for all practical purposes to the results we
would get if we used a truly instantaneous process.
We could use some other interest rate than 100% in the calculations just presented. In fact we could use any interest rate at all. Similarly we could use any time period over which to calculate the cumulative value of our investment. Our rate will be in terms of some time unit. The n or number of compoundings now refers to that particular time unit. We did not have to stop our compounding process at one time unit. Given very small compounding steps, we could stop at almost any point before one time interval or go for any period of time afterwards. Similarly, we did not have to start out with only one unit of money on which we were accumulating interest. The general mathematical representation for the amount of money accumulated by leaving an original amount {M(time 0)} to accumulate interest over an interval "t" time units long at a rate "r" of interest when interest is compounded "n" times over one time unit is:
....................................Equation(5)
Here "t" is the number of time units over which we will
be accumulating interest. "n" is the number of times
the interest will be compounded during each one of those time
periods. "r" is the interest rate. 
is the amount of money that we start with. Go through the logic
for this equation like we did for the 100% rate and 1 time unit
period case to be sure you understand where this equation comes
from.
Let us state again what we mean by a closed form solution of a dynamical system and contrast that to a simulation of a dynamical system. If we have a closed form solution and we know the values of the entities in a system at one point in time, we can plug the value of any time interval into our closed form solution and calculate the values of the entities of the system after that time interval. With simulations, on the other hand, in order to determine the values of the entities of the system at any later point in time, the simulation must take tiny steps from the starting time until it reaches the time of interest. With equation (2) above, all we need to calculate the value of M at time t is the value of m at time zero, the interest rate r, the value of n we will use, and a good calculator. If we had n approaching infinity, it would seem as if we were approaching the closed form solution for differential equations because with n approaching infinity, the time step approaches zero. The question is, however, would your calculator just be making all the simulation steps?
Before we find a more compact closed solution that our calculator might use, lets look at the effect of increasing the number of compoundings per year (n) when our time interval is three years and our interest rate is 50% per year and our starting amount is $100. The following table 2 presents the calculations.
| Number of compoundings per year | Formula for accumulated value after three years | Accumulated value after three years |
| 1 | 100(1+0.5/1)^1*3 | 337.5 |
| 2 | 100(1+0.5/2)^2*3 | 381.5 |
| 4 | 100(1+0.5/4)^4*3 | 411.0 |
| 8 | 100(1+0.5/8)^8*3 | 428.4 |
| 16 | 100(1+0.5/16)^16*3 | 438.0 |
| 32 | 100(1+0.5/32)^32*3 | 443.0 |
| 64 | 100(1+0.5/64)^64*3 | 445.6 |
| 128 | 100(1+0.5/128)^128*3 | 446.9 |
| 256 | 100(1+0.5/256)^256*3 | 447.5 |
| 512 | 100(1+0.5/512)^512*3 | 447.8 |
| 1024 | 100(1+0.5/1024)^1024*3 | 448.0 |
| 2048 | 100(1+0.5/2048)^2048*3 | 448.1 |
| 4096 | 100(1+0.5/4096)^4096*3 | 448.1 |
| 8192 | 100(1+0.5/8192)^8192*3 | 448.1 |
| 16384 | 100(1+0.5/16384)^16384*3 | 448.2 |
Describe step by step how you could get
Stella® to make these calculations. Do the calculation for
500 both on a spreadsheet or calculator and using STELLA®.
Note that when the compoundings get to very high numbers, the
dt required in the Stella model becomes too small. In that case,
you are going to have to change the units of the rate parameter
to a smaller unit. For example, instead of using 0.5 per year,
you may have to use 0.0416666666666667 per month.
Let us reconsider how each of the calculations in table 1 was
made. We will take "n"=4 as an example. At "n"=4,
there are 12 compoundings, one every three months, over the 3
year period. At each compounding we take 0.5/4 or 12.5% of the
value at the beginning of the period and add it to the original
value. That is like taking the original value and multiplying
it by 1.125. Then taking that value and multiplying it by 1.125.
Then ...etc. for a total of 12 times. Thus we have 100*1.125^12.
We can plug that into a calculator.
Instead of working out the formula so that each additional multiplication
adds one to the value of the exponent in formula one, lets work
out the calculations like a simulation does, step by step. The
following table was generated on a spreadsheet. Make sure that
you can relate these calculations to those in a Stella® table
with each of the settings for "Report" and "Report
flow values" in the table control window.
In table 3 we are simulating an exponential growth process just
as Stella® would using the Euler method. In an exponential
growth process, the increase in size over a very short interval
is proportional to the size of the population. The following graphs
are plots of this growth process first just over three years and
using exactly the data in table 2. Then the same process is extended
over a 30 year period. 
Note that the data in the first 3 years of the second figure is
the same as the data from the first 3 years of the first figure.
When you extend the time period on an exponential growth process,
it looks pretty explosive. The point of showing these two graphs
is to make the point that every exponential process is ultimately
explosive, even if over a short time span it might not seem so.
Because population growth was an exponential process, some people
projected that the earth's population would explode like a population
bomb. If we extend the current population growth processes, we
could see that it wouldn't be too long before the mass of humans
on earth would be greater than the mass of the earth itself. Of
course exponential growth couldn't possibly continue to that point.
You need to be sure that you can see how this process in table
3 relates to equation (5) above.
If "r" and "t" and 
are other than "1" the limit reached as 
is 
.
To see this, define "m" as 
.
Then 
.
As n tends to infinity, m tends to infinity. The limit of 
as m approaches infinity is "e" by definition
and therefore
.....................................Equation
(6)
Any quantity to the zero power is one so that at t = 0, there
is M(time 0) in the account.
Note also from the above examples that
when rt is very small.
That is to say that the simulated value across a small discrete step is approximately the same as the value from the closed form solution of the differential equations. We will call the small discrete steps "dt". When "n" is large, dt is small. When we our doing simulations using STELLA™, we can get the simulation closer to a closed form solution of the differential equations by making "dt" smaller.
Note that with equation (6), in order to get our accumulated value at any point in time, we do not need to go through the kind of calculations we did in table two as long as we have a calculator with a natural logarithm function on it. It turns out, however, that to make the calculation of the base of the natural logarithm, some approximating model is still used by the calculator. The calculations in table 3 correspond to what we will call simulation of a dynamic process using a difference equation approximation to the differential equations. When we have a closed form solution like equation (6), we can plug the time and the rate and the initial values into the formula and get our accumulated value at any point in time directly in a way that more directly corresponds to differential equations rather than their approximation by difference equations.
In the above example, the population was money. But it could have
been any population of individuals whose size was being changed
at a constant rate proportional to its current size at any point
in time. Another generalization is that instead of growing, one
could be losing money or population at a steady rate. That is
to say the rate could have a negative value. We are now going
to generalize the above conclusions about population size given
exponential processes to negative growth rates. In the Stella™
model of a fixed rate of disease in a population which we presented
earlier, we had a negative compounding of the number of well individuals.
By treating the rate of disease development as a negative growth
in a healthy population, we are going to derive the relationship
between disease rates and disease risks.
Say that the rate of a disease was 1/100 per month (1 case per
100 individuals per month or 1 case per 100 person months) and
that once one gets the disease they are never healthy again. Then
after one month the size of the healthy population would be smaller
than when it started. A first approximation as to how much smaller
W would be is
. The remaining well population
would be 
. This calculation assumes that
the rate of per month can be treated like a risk over one month.
We will have more to say later as to when this will be a reasonable
assumption.
Note that when we subtract a proportion from a population rather
than adding a proportion to a population, the general formula
for the size of the population after a number of compoundings
is very similar to equation (5). It is
...........................Equation (7)
As we take the limit as n tends to infinity using the same approach
as with equation (2) we get
..................................Equation
(8)
Equations(7) and (8) are the same as equations (5) and (6) but
just with negative growth rates rather than positive growth rates.
The risk of a disease results from a rate of disease acting on
a population over a period of time. We saw that the proportion
of a population remaining healthy when a disease rate was acting
on it over a time period t was 
. The proportion
of the population getting the disease, that is to say the risk
of the disease is one minus that. In other words:
More precisely we might say
This is a fundamental relationship upon which a lot of epidemiologic
theory is built. It embodies models of risk and models of rares.
It assumes that the rate stays constant over the entire time period.
For some approximation purposes we go ahead and make this assumption
about the constancy of the rate. In most real dynamic systems,
the rate will be a function of other variables in the system and
will vary over time.
We intend to show now that risks are always numerically less than
the rates which generated them when the time unit in the denominator
of the rate is the period across which the risks are calculated.
That is to say that a risk across time "t" is numerically
less than a rate per person per time "t". We will demonstrate,
however, that risks of 0.05 or less can for most practical purposes
be translated directly into constant rates because they are only
insignificantly smaller than the rates. This near equivalence
of risks and rates under the condition of small risks may be one
reason why risks and rates are used interchangeably by many authors.
The situations where risks and rates are not interchangeable,
however, are frequent enough so that we have to make a clear distinction
between risks and rates.
The reason why the numerical value of a risk is less than that
of a rate if both use the same time period is that a rate acts
instantaneously so that as soon as one person becomes ill, the
rate is in fact acting on a smaller population. We think of a
rate per time period as acting continuously. That is to say that
in any short interval where we could make an observation, we would
expect the rate to be acting to generate cases and removing healthy
individuals from the population base on which the rate acts. This
of course would only be the case if the population is very large.
When we model populations using STELLA, this assumption about
the continuous action of rates is thus almost the same as assuming
that a population is very large.
Given a very large population, we could expect the size of the
population to change in a relatively smooth fashion instead of
a jerky fashion as would be the case if the population size were
10. If the population were that small, when one individual became
ill we would have an abrupt jump of the population size to 90%
of its previous value. Given a large population we could calculate
the real size of the population on which the rate is acting at
an increasingly large number of interval points over the month.
If our disease rate is 0.01 per month, in the limit then the fraction
of the original population remaining after a month would be
M(time t) = e-rtM(time
0) = e-0.01M(time 0)
= 0.99004983M(time
0).
The situation where we used the rate as if it were a risk, (that
is the approximation that assumes the population on which the
one tenth per month rate acts on stays constant over the month)
differs from the case where the rate is used properly by only
an insignificant amount (0.99 vs 0.99004983).
That is to say at these values of "r" and "t",
there is not much difference in the value of a rate and a risk.
Even if we have a rate of 0.1 per time period instead of 0.01,
there would not be all that much difference between the numerical
value of a rate and the risk over which that rate acts during
one time period. e-.1 = .09484 which is not
so far from 0.1. As you accumulate this kind of error over more
steps, however, it would begin making a difference. To get a feel
for how constant rates act to generate cumulative risks, we present
below the step by step cumulative process of a rate of 0.01 per
time period and a rate of 0.1 per time period acting over 25 time
periods. We also present the exact closed form solution for the
risk acting over the 25 time periods using the natural logarithm
exponentiation formulas presented above.
| Months | Remaining populations | Remaining populations | |||
| "Rate" =.1 in a simulation where the rate is used as if it were a risk. | Rate = .1 closed form | Rate = .01 in a simulation where the rate is used as if it were a risk. | Rate = .01 closed form | ||
| 0 | 1 | 1 | 1 | 1 | |
| 1 | 0.9 | 0.904837418 | 0.99 | 0.990049834 | |
| 2 | 0.81 | 0.818730753 | 0.9801 | 0.980198673 | |
| 3 | 0.729 | 0.740818221 | 0.970299 | 0.970445534 | |
| 4 | 0.6561 | 0.670320046 | 0.96059601 | 0.960789439 | |
| 5 | 0.59049 | 0.60653066 | 0.95099005 | 0.951229425 | |
| 6 | 0.531441 | 0.548811636 | 0.941480149 | 0.941764534 | |
| 7 | 0.4782969 | 0.496585304 | 0.932065348 | 0.93239382 | |
| 8 | 0.43046721 | 0.449328964 | 0.922744694 | 0.923116346 | |
| 9 | 0.387420489 | 0.40656966 | 0.913517247 | 0.913931185 | |
| 10 | 0.34867844 | 0.367879441 | 0.904382075 | 0.904837418 | |
| 11 | 0.313810596 | 0.332871084 | 0.895338254 | 0.895834135 | |
| 12 | 0.282429536 | 0.301194212 | 0.886384872 | 0.886920437 | |
| 13 | 0.254186583 | 0.272531793 | 0.877521023 | 0.878095431 | |
| 14 | 0.228767925 | 0.246596964 | 0.868745813 | 0.869358235 | |
| 15 | 0.205891132 | 0.22313016 | 0.860058355 | 0.860707976 | |
| 16 | 0.185302019 | 0.201896518 | 0.851457771 | 0.852143789 | |
| 17 | 0.166771817 | 0.182683524 | 0.842943193 | 0.843664817 | |
| 18 | 0.150094635 | 0.165298888 | 0.834513761 | 0.835270211 | |
| 19 | 0.135085172 | 0.149568619 | 0.826168624 | 0.826959134 | |
| 20 | 0.121576655 | 0.135335283 | 0.817906938 | 0.818730753 | |
| 21 | 0.109418989 | 0.122456428 | 0.809727868 | 0.810584246 | |
| 22 | 0.09847709 | 0.110803158 | 0.80163059 | 0.802518798 | |
| 23 | 0.088629381 | 0.100258844 | 0.793614284 | 0.794533603 | |
| 24 | 0.079766443 | 0.090717953 | 0.785678141 | 0.786627861 | |
| 25 | 0.071789799 | 0.082084999 | 0.777821359 | 0.778800783 | |
After 25 time periods of using a rate of 0.1 as if it were a risk,
the simulated value of 0.071789799
is about 12% less than the more exact closed form calculation
which gives a value of 0.082084999.
After 25 time periods of using a rate of 0.01 as if it were a
risk, however, the difference from the exact closed form solution
is less than 0.1% (0.777821359 vs.
0.778800783).
Find the risk of a disease over a 12 month period of a disease that is occurring at the rate of 0.04 per person per month. Do the calculation first using natural logarithms and then using the approximation where you make 12 monthly risk calculations. Next do the calculation using 48 0.25 month steps over the one year period. Use both a spread sheet and Stella™ to make the simulation calculations. Integrate both into your word processing document.
A general rule of thumb is that a rate can be used as if it were
a risk if the value of the risk is 0.05 or less over the interval
considered. That same rule of thumb does not apply, however, to
simulations where the errors generated across one step can be
compounded in the next step. We will see that in doing simulations,
we need to check to see that the end result isn't altered by the
step size we take.
There are various situations in which we may want to use a risk
as a rate or vice versa. In this class, we will be wanting to
enter rates into simulations to examine dynamic consequences.
The field data we have relevant to rates will sometimes be risk
calculations over some defined period. That is to say, the field
studies will not have determined the actual person-time in the
denominator for the rate. They will just have used the number
of people at risk at the beginning of the period as their denominator.
The true rate that generated the risk is always greater than the
observed risk because its denominator is really smaller than the
denominator used in the risk.
When can we treat those risk calculations as rates to be entered
into our calculations? When for all practical purposes the numerical
differences between the risk over a unit of time is not meaningfully
different from the numerical value of the rate. That will only
be the case when the risk is small. How small? When we are running
simulations, we will cut our step times until cutting the step
times makes no further difference. At that point the risk is small
enough. In presenting reports, a good rule of thumb is that for
all practical purposes a risk of 0.05 or less can be treated as
a rate per the time unit over which the risk was calculated. What
that means in doing simulations in STELLA™ is that you do
not want the value of any of your stocks to change more than 5%
across any DT.
Sometimes we will be given a rate where the person-time in the
denominator was accurately determined. We may want to use that
rate to express a risk to an individual. The actual risk will
always be numerically lower than the corresponding rate. When
can we use the rate as a risk? Again the rule of thumb of a risk
< 0.05 can be accepted for most practical purposes. Sometimes
we want to express a risk over a longer time period than the period
in which the rate is expressed. For example we might want to express
a risk across ten years when we know that the rate per year is
0.02. Assuming constant rates, when can you just multiply the
rate times the time to get the risk? Another way of asking this
question is, "When can we for practical purposes disregard
the fact that as a rate is extended over time it is acting on
fewer and fewer unaffected individuals and just treat the denominator
as if it is staying constant? Again the rule of thumb that when
your risk estimate comes out to be less than 5% this might be
OK. In other words, when rate times time is less than 0.05.
Note that in the kind of simulation we have constructed, the time when a case entered the well population has nothing to do with the time when it flows out of the well population. The reservoir stock of a STELLA™ model assumes that all cases in a stock become thoroughly mixed and that the flow out of the stock has nothing to do with the timing of the flow into the stock. This independence from the timing of previous events is sometimes called the Markov assumption. This assumption is very convenient for computer simulations. The computer does not have to keep track of what has happened in the past when this assumption is made. This is a characteristic of ordinary first order difference equations. That is the kind of equations that STELLA™ uses.
Whenever the flow out of a state depends on what has happened
to a system in the past, there is a dramatic increase in the amount
of computer work that must be done. For each flow in the model,
the previous states of the system have to be reviewed and calculations
made on the basis of these. Just storing the previous states may
overwhelm the usual desk top computer. Doing all the calculations
based on the previous states can slow down simulations so that
they might take days. Thus it is to our advantage to learn to
work with the Markov assumption.
Epidemiology 802 Course Home Page.