One of the current social issues receiving a fair amount of publicity in the media is the question of the age distribution of the human population in the United States. The usual context for this problem is the question of the future solvency of the social security system, but there are many other important questions dealing with future populations and the age structure of that population. (The need for new elementary, middle and high schools, for example.)
In this section of the tutorial you will build a simple model which forecasts an age structured population many years into the future. The assumptions are based on the premise that certain demographic parameters remain constant over the period of the time being studied. As the tutorial progresses, I will develop an example of this type of model called the Leslie Matrix model while encouraging you to develop a similar model on your own.
The main idea of these models is to split a population into several distinct age classes and model how individuals are created and subsequently move from one age class to the next. In the Leslie model every member of one age class either dies or moves into the next age class at every iteration of the model. In the model you will be constructing individuals will experience one of three possibilities: death, stasis (remaining in the same age class) and maturation (moving on to the next age class).
Focusing for awhile on the Leslie model, there are two parameters associated with each age class. The percent survival rate for the age class represents the proportion of the individuals in an age class which survive to move into the next age class during the next iteration of the model. The second parameter is the fecundity parameter which represents how many offspring (age class zero) are produced, on average, by each member of the age class. (Some argue convincingly that fecundity is most closely related to the number of females in the population and thus create Leslie models which count only the females in a population.)
To continue our development of the Leslie Model we require an
introduction of some mathematical notation. Suppose P(i)
represents the population of age class i. If s(i)
represents the survival rate of age class i, then the
population in age class i+1 in the next iteration can be
computed using the current population in age class i through
the formula
next iteration's P(i+1) = s(i)*P(i).
If we let N be the last age class in our model (all
individuals in age class N die before the next iteration)
then the formula above works for values of i from 0 up to
N-1, inclusive. This still leaves us without a population
for the next iteration's youngest age class, age class 0.
The fecundity parameters are used to determine how many newborns are
created during the next iteration of the model. If we let
f(i) be the fecundity (average birth rate) per individual in
age class i, then the total number of newborns will be:
f(0)*P(0) + f(1)*P(1) + f(2)*P(2) + ... + f(N)*P(N)
If P is a column vector of population by age class, then it is possible to compute the P vector for the next iteration through a simple matrix (non-componentwise) multiplication. In particular, if A is the (N+1)x(N+1) matrix whose first row is the fecundity rates f(i) and whose diagonal below the main diagonal is made up of the survival rates s(i), then the next iteration's population column vector can be found by performing the matrix multiplication A*P.
As an example, suppose we have a population of critters who never grow past five years old. Suppose, further, that 80% of the critters survive each year to move onto the next age class, with the exception of the five year old critters which always die. If newborns and yearlings have no offspring, and if middle aged critters (2, 3 and 4 year olds) have, on average, 0.35 newborns each year, while the oldest critters (5 year olds) have only 0.1 newborns on average, how can we use current populations to forecast future populations and the age structure of that population?
In this model the parameters are given in the table below:
Age Class Survival Rate Birthrate
i s(i) f(i)
0 0.8 0.00
1 0.8 0.00
2 0.8 0.35
3 0.8 0.35
4 0.8 0.35
5 --- 0.10
To create the matrix A mentioned above for this particular set of parameters, the MATLAB command below would do the trick. (Note the elipses at the end of some lines. This is how we enter a MATLAB command which takes more than one line to complete.)
A = [ 0.00 0.00 0.35 0.35 0.35 0.10; 0.8 0 0 0 0 0; ... 0 0.8 0 0 0 0; 0 0 0.8 0 0 0; 0 0 0 0.8 0 0; ... 0 0 0 0 0.8 0]which creates the matrix A shown below:
A = 0.0000 0.0000 0.3500 0.3500 0.3500 0.1000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.8000 0.0000
To forecast the future age structured population of our colony of critters we need only determine the current age structured population, represent this population as a column vector P and then move forward one year at a time in our model by performing the multiplication A*P.
For example, if our initial population consists of sixty critters evenly spread throughout the six age classes, our initial P vector would be:
P = [10;10;10;10;10;10]Note that the semicolons cause P to be a column vector. We could have achieved the same end by separating the 10's with spaces or commas, creating a row vector, and then using the matrix transform operator (an apostrophe).
To determine the population distribution next year we multiply A and P giving the vector
11.5000 8.0000 8.0000 8.0000 8.0000 8.0000If we want the population two years in the future there are a number of ways to proceed. First, we can take the vector above and multiply it by A to get the population for the next iteration. Proceeding in this fashion we get a series of vectors, each of which gives the population for the succeeding iteration.
If we only want populations for specific future years a second way may make more sense. This method is based on noting that for each year into the future we go, we arrive at the population by multiplying by the matrix A. Thus, if we want to go ten years into the future we get there by multiplying by A ten times. We can simplify this process by using the (noncomponentwise) matrix exponentiation operator carat (^).
To move ten years into the future we start with the current population (held in the column vector P) and multiply it by A ten times. This could be accomplished with the MATLAB command
A*A*A*A*A*A*A*A*A*A*PThere is an easier way however. Note that in the computation above, we multiply the matrix A by itself 10 times. (This is traditional matrix multiplication, not componentwise multiplication.) We can use the exponention operator (^) as a way to indicate repeated matrix multiplications without writing them all out. Using matrix exponentiation, the MATLAB command given above can be recast as:
A^10*PIt is important to note that the carat (^) without a leading dot indicates repeated traditional matrix multiplication. This contrasts with the dot-carat (.^) version of the exponentiation operator which would work in a componentwise fashion, raising each entry in a matrix to the indicated power.
In our example, the command A^10*P yields the result:
3.0288 2.7583 2.4885 2.3087 2.1601 1.9545a result which leads us to be very concerned about the long term viability of this population.
ones(1,6)*(A^10*P)
To investigate this question we must update the entries in the A matrix to reflect the hypothetical introduction of the new drug. Only three entries inthe matrix A need be changed, so rather than recreate A from scratch we will make the three individual changes.
The first change is that the birth rate for two year olds must be changed to 2.3. Since the birth rate for two year olds is found in the first row, third column of the matrix A, this change can be effected by the MATLAB command:
A(1,3)=2.3;
Make the remaining two changes in A, changes which reflected the lowered survival rates for age classes 3 and 4.
Begin with the same initial population (all 10's) as before. Does the population seem healthier (whatever that means) under the drug regime? How does the drug effect the age distribution of the population?
If fecundity rates remain unchanged, what survival rates are needed for our critter population to remain stable? What proportion of a stable population falls into each age class?
It is fun to see how populations evolve over time. Toward this end it is possible to create a three dimensional graph showing the evolution of a Leslie population over time by creating a matrix whose columns are "snapshots" of the population at various times under study.
Suppose P is the initial population (all 10's in our example). We will use the matrix TimePop to be the matrix whose columns will contain the population at given points in time. Since the first population we are interested in is the initial population, we will set the first column of TimePop to the initial population using the MATLAB command
TimePop=P;We will use a second matrix called PP to hold the most recent population computed. Thus, PP should start off holding the initial population
PP=P;The task now facing us is to generate the population vector for the next time period, add this column vector to the matrix TimePop and then do the same again, moving ahead yet another time step. This can be achieved by executing the following MATLAB commands repeatedly:
PP = A * PP; TimePop = [TimePop PP];The first command moves the population vector PP forward one time step. The second command adds a new column to the TimePop matrix. This new column contains the most recently computed population.
Using the up-arrow key, repeat the line above 15 to 20 times, creating a TimePop matrix with many population snapshots.
The evolution of the critter population over time can be viewed by creating a three dimensional surface plot of the array TimePop. To create this surface plot enter the following MATLAB commands and note the plot which is generated. Pay special attention to the creation of labels and the colorbar.
surf(TimePop)
view(0,90)
colormap(jet)
colorbar
xlabel('Year')
ylabel('Age Class')
The colon (:) notation we discussed earlier is a convenient way to create large arrays where the entries are evenly spaced. A second use of the notation is to extract parts of existing matrices. In particular, a single colon is a reference to every row or every column in a matrix, depending on its position.
For example, the notation A(2,:) refers to every column of the second row of the matrix A, or, said another way, A(2,:) is the entire second row of A.
Similarly, B(:,5) refers to the entire fifth column of the matrix B.
plot(TimePop(3,:))
Using calculators it is difficult to efficiently deal with Leslie models consisting of more than about ten age classes. Using MATLAB, however, it is relatively easy to deal with populations with dozens of age classes. Human populations are generally modeled with 101 age classes, allowing for ages between 0 and 100. It is also much easier to move far into the future, spotting long term trends.
In this section we discuss modifications to the Leslie Population Model which describe a slightly different situation than that modeled by Leslie. This section introduces no new MATLAB skills and is present for those who want to undertake the challenge of creating a new model. I will describe the mathematics of the new model, but will leave it entirely up to you to implement the model in MATLAB. A knowledge of traditional matrix multplication, the non-componentwise variety, will be helpful.
In the Leslie model only two things can happen to individuals in an age class; they can die or they can move onto the next age class. In the model presented here we allow for a third option, remaining alive and in the same age class.
Given that it is possible to remain in the same class for several time steps, it is probably best to leave off the word age when describing these classes.
As far as I know, this version of the Leslie model has no well-known name, so we get to name the model for the purposes of this tutorial. Feeling singularly uncreative at the moment, I choose to call it Model X.
Model X requires three sets of parameters. Denote the fecundity of class i with the symbol f(i), the proportion of class i moving on to the next class at each iteration with the symbol s(i) and the proportion of each class which survives but remains behind with the symbol r(i).
Proceeding in a fashion reminiscent of the Leslie model development, if we create a model with classes numbered from zero to N, and let P(i) denote the number of individuals in class i, then we can derive the following equations which give the values of the various P(i) for the next time step in terms of the P(i) values for the current time step.
Considering, first, the values of P(i) in the next time step
for classes numbered from one to N we see:
Next period's P(i+1) = s(i)*P(i) + r(i+1)*P(i+1)
This equation tells us that next year's class i+1 is made up
of those individuals from this year's class i who move up a
class as well as those individuals from this year's class i+1
who survive but do not advance in class. (a.k.a. flunking?) This
equation holds for values of i from zero to N-1
inclusive.
A careful look at this situation reveals that we have not yet specified the population of class zero in next year's population. Recall that members of class zero are either true newborns or previous members of class zero who do not advance for some reason.
The mathematical expression which gives the number of class zero
residents in next year's population is given by:
Next period's P(0)=
(r(0)+f(0))*P(0) + f(1)*P(1) + f(2)*P(2) + ... + f(N)*P(N)
Pay special attention to the first term of this equation. It does not
follow the pattern of the other terms.
Model X has a structure similar to the Leslie model. In fact, if P is a column vector whose N+1 entries are the populations of each class at some point in time, then there is an (N+1)x(N+1) matrix A so that the population next year can be determined by calculating the traditional matrix product
A*P