This post continues my tutorial revisiting the famous “Artificial Anasazi” model of Axtell et al. (2002), ported to NetLogo by Marco Janssen in 2008 (see his paper). I have rewritten the model for explanatory purposes (though very much inspired by Janssen’s work). Part 1 of the tutorial looked at the patch dynamics, and we now turn to modelling the people.
The model studies the population of Anasazi farmers living in Long house Valley (near Black Mesa in Arizona) between 800 and 1350 AD. The image below shows the NetLogo view of this valley, as modelled for 1000 AD, with each person representing a household. For purposes of simplicity, I am not modelling the availability of drinking water and the locations of dwellings. The locations of the people on the diagram below therefore represent farm locations. However, because I am not modelling water availability and dwelling locations, the spatial distribution of people should not be taken too seriously – I am only hoping to accurately model population numbers.
Anasazi society was matriarchal (households were headed by women) and matrilocal (daughters lived near their mothers). Since archaeology (counting and dating house ruins) gives us household numbers (rather than people numbers), the model operates at a household level. However, food consumption calculations are based on five people per household. Young women are assumed to marry at age 16, so that the possible history of a household is as follows:
Household age | Matriarch age | First possible daughter | Last possible daughter |
---|---|---|---|
0 (begins) | 16 | – | – |
1 | 17 | 0 (born) | – |
17 | 33 | 16 (marries) | – |
18 | 34 | – | – |
19 | 35 | – | – |
20 | 36 | – | 0 (born) |
36 | 52 | – | 16 (marries) |
37 | 53 | – | – |
38 | 54 | – | – |
The earliest age that a household can reproduce (spin off a daughter household) is 17 (at which time the mother is aged 33). Fertility is assumed to continue to age 36 (based on the discussion in Axtell et al.), which means that the latest possible daughter household is produced at a household age of 36. Matriarchs are assumed to die after age 54, so that the oldest possible household is 38 years (which is perhaps a little old for a society of this kind).
The Anasazi were able to store corn (according to Axtell et al., for 2 years), which gave them the ability to ride out short droughts. Daughter households were gifted some of their mother’s corn (about one third), to get them started in life. This gives us the following constants for the model:
to setup
...
set start-year 800 ;; initial year
set food-requirement (5 * 160) ;; a person needs 160 kg of corn each year, while a typical household consist of 5 persons
set fertility-age 17 ;; the minimum age of agents (households) that can reproduce (spin off daughter households)
set fertility-end-age 36 ;; the maximum age of agents (households) that can reproduce
set death-age 38 ;; the age after which agents (households) die (dissolve when the matriarch dies)
set corn-gift-ratio 0.33 ;; each new household gets this fraction of the corn storage of the parent household
set corn-stock-years 2 ;; corn can be stored for this many years
...
end
Households must make certain decisions. Initially, and in times of famine, they must choose a (new) farm. We introduce a reporter (function) to return an estimate of the productivity of a farm (patch). This reporter will be used twice, so it is good practice to give the calculation its own name. The result of the reporter is in fact just be the base yield of a patch, but we could expand it to incorporate historical information and/or decision errors. Giving the calculation its own name makes it easier to include such extensions down the track.
to-report estimated-farm-yield
report base-yield
end
We use this estimate to generate an agent set of potential farms – those that produce enough food for a household and are not being farmed by anybody else:
to-report potential-farms
report (patches with [ estimated-farm-yield >= food-requirement and not being-farmed? ])
end
The following reporter picks out the best farm from an agent set. If there are no suitable farms, the household is assumed to leave the valley, which is equivalent (from a model point of view) to dying. But what is the “best farm”? The matrilocal structure means that agents try not to move very far (i.e. they select a low distance from their existing location), but we also introduce a slight bias towards a high estimated farm yield (minimising 1000 / estimated-farm-yield
). The NetLogo built in min-one-of
operator does the hard work of finding the patch with the lowest combined score. Incidentally, for distances to be calculated correctly, the NetLogo world must be set to not “wrap” either horizontally or vertically.
to find-best-farm [ farms ]
ifelse (count farms = 0)
[ die ]
[ let existing-farm patch-here
let best-farm (min-one-of farms [ distance existing-farm + 1000 / estimated-farm-yield ])
ask best-farm [ set being-farmed? true ]
move-to best-farm ]
end
According to the archaeological record, the population of the valley began with 14 households. We can initialise them with the code below (which assumes that households are what NetLogo calls a “breed”). We use a utility procedure initialise-household
to set various attributes. This procedure takes a starting age as a parameter. Starting ages run from 0 to 28, but we will use the same procedure later for new-born households, which will always start with an age of zero. We also give each household an initial list of corn storage stocks, divided into ages of 2, 1, and 0 years. The oldest corn (at the front of the list) will be eaten first, and corn older than 2 years will be thrown away. The NetLogo built in n-values
operator generates a list of the required length, calling random-float
three times.
to setup
...
set-default-shape households "person"
calculate-patch-yields ;; needed to set up decision-making
create-households 14 [
initialise-household (random 29)
set corn-storage (n-values (corn-stock-years + 1) [ 600 + random-float 400 ])
find-best-farm potential-farms
]
...
end
to initialise-household [ start-age ]
set age start-age
set harvest 0
set size 4
set color black
end
Closely related code allows a household to reproduce. For efficiency reasons, we pass the set of potential farms as a parameter. The NetLogo built in hatch
operator allows an agent to create other agents. We also use the NetLogo built in map
operator twice – once to generate a list of gift corn, with each entry being one third of the parent household’s list; and once to subtract the entries in the gift-corn list from the corresponding parent household’s list (as a general rule, if a list is being transformed to another list of the same length, the map
operator usually provides the most elegant NetLogo solution).
to do-reproduce [ farms ]
let gift-corn (map [ ? * corn-gift-ratio ] corn-storage)
set corn-storage (map [ ?1 - ?2 ] corn-storage gift-corn)
hatch 1 [
initialise-household 0
set corn-storage gift-corn
find-best-farm farms
]
end
We can now provide the code for a simulation step, although annual-household-activities
and plot-interesting-data
still need to be defined. There are two ask
loops. In the first loop, agents set off in search of a new farm if they predict that their corn stocks (plus their expected future harvest) will not meet their annual food requirements (plus a 10% margin). In the second loop, households of fertile age reproduce with a probability given by the fertility
slider. A local variable stores the result of the potential-farms
calculation, to avoid recomputing it unnecessarily.
to go
calculate-patch-yields
annual-household-activities
ask households [ ;; agents with insufficient food for next year try to find a new farm
if (food-estimate < 1.1 * food-requirement) [
ask patch-here [ set being-farmed? false ]
find-best-farm potential-farms
]
]
let farms potential-farms ;; avoid recomputing this unnecessarily
ask households [
if (age >= fertility-age and age <= fertility-end-age and random-float 1 < fertility and count farms > 0) [
do-reproduce farms
set farms potential-farms
]
]
plot-interesting-data
if (year = 1350 or count households = 0) [ stop ]
set year year + 1
tick
end
to-report food-estimate ;; estimate the amount of food available for next year, based on current stocks of corn, and an estimate of the future harvest
let future-estimate harvest ;; predict next year's harvest to be the same as this year's harvest
report (future-estimate + sum (but-first corn-storage)) ;; ignore the corn on the front of the list, which is now too old to eat
end
The annual-household-activities
procedure increments the age of households, has them eat their annual food quota, and kills them off if they are too old or have insufficient food. Points to note include:
- We set up a global variable
harvest-list
to be a list of this year’s harvests (this will later be plotted as a histogram). - We supplement the spatial variability in farm quality with a temporal variability in annual harvests, re-using our
apply-variability
reporter. - We maintain the household’s 3-element list of corn by dropping the oldest corn at the front, and adding the new harvest to the back.
- We introduce a household attribute
unsatisfied-hunger
to process the food that the household needs. - We use a NetLogo trick to eat corn starting from the front of the list, using the
eat-corn
reporter (below).
to annual-household-activities
set harvest-list [ ] ;; a list of this year's harvests (1)
ask households [
set age age + 1
if (age > death-age) [
ask patch-here [ set being-farmed? false ]
die
]
set harvest (apply-variability harvest-variability) * ([ base-yield ] of patch-here) ;; this year's harvest (2)
set harvest-list (fput harvest harvest-list) ;; (1)
set corn-storage (lput harvest (but-first corn-storage)) ;; oldest corn is at the front of the list (3)
set unsatisfied-hunger food-requirement ;; a household attribute (4)
set corn-storage (map [ eat-corn ? ] corn-storage) ;; eat corn starting from the front of the list (5)
if (unsatisfied-hunger > 0) [
ask patch-here [ set being-farmed? false ]
die
]
]
end
to setup
...
set harvest-list [ ]
...
end
In transforming a list to another of the same length, using map
is usually a good idea. Here it processes the list from the front, applying eat-corn
to each element. The unsatisfied-hunger
attribute links together the different calls to eat-corn
, so that they either (a) subtract the unsatisfied food requirement from a list element, or (b) consume all of a list element and update unsatisfied-hunger
for the next iteration. The eat-corn
reporter thus has two simple cases:
to-report eat-corn [ existing-amount ]
ifelse (existing-amount >= unsatisfied-hunger)
[ let remaining-amount (existing-amount - unsatisfied-hunger)
set unsatisfied-hunger 0
report remaining-amount ]
[ set unsatisfied-hunger (unsatisfied-hunger - existing-amount)
report 0 ]
end
A diagram showing three example calls to eat-corn
may make this clearer. It can be seen that the oldest corn at the front of the list is being eaten first:
The following code handles plotting and statistics. We use plotxy
so that our x-coordinates can start at 800. We plot histograms of household ages and our previously calculated harvest-list
(we don’t use the obvious [ harvest ] of households
because that would exclude households that starved to death, and would also include new daughter households that have no harvest yet). Unfortunately, histograms in NetLogo require explicit control of the x axis, which adds some complexity. We also maintain a list of population levels in population-sequence
, and calculate how well this fits the historical sequence (excluding the proportion of the historical sequence where the population is zero). We measure fitness using the root-mean-square of differences between the two sequences (the aquare root of the mean of the squares of the differences). This is equivalent to using what the literature calls the L2 norm (which penalises large errors), but is expressed in units of numbers of households, which is easier to understand. Values of the root-mean-square error tend to be around 30, which means that “typically” the model differs from the archaeological reality by about 30 households.
to plot-interesting-data
set-current-plot "Population (households)"
set-current-plot-pen "Potential"
plotxy year (count patches with [ mean historical-base-yield >= food-requirement ])
set-current-plot-pen "Historical"
plotxy year (item (year - 800) historical-sequence) ;; index list starting at 0
set-current-plot-pen "Model"
plotxy year (count households)
let age-list ([ age ] of households)
set-current-plot "Household ages"
let xmax (1 + max (fput 0 age-list))
set-histogram-num-bars xmax
if (plot-x-max < xmax) [ set-plot-x-range 0 xmax ] ;; histograms require explicit x-max control
histogram age-list
set-current-plot "Household harvests"
let block-size 100
set xmax (1 + ceiling (max (fput 0 harvest-list) / block-size))
if (plot-x-max < xmax * block-size) [ set-plot-x-range 0 (xmax * block-size) ]
set-histogram-num-bars (plot-x-max / block-size)
histogram harvest-list
if (length population-sequence < desired-fit-length) [ set population-sequence (lput (count households) population-sequence) ]
let cut-historical-sequence (n-values (length population-sequence) [ item ? historical-sequence ])
set fit-badness (sqrt (mean (map [ (?1 - ?2) * (?1 - ?2) ] population-sequence cut-historical-sequence)))
end
to setup
...
set population-sequence [ ]
set desired-fit-length (length filter [ ? > 0 ] historical-sequence)
...
end
A typical plot is shown below (simulated population levels are in red). The root-mean-square error here is 25.963, which is quite good. Our decision to introduce the usable-farm-fraction
slider is partially responsible for this, since it limits excessively high populations without unduly restricting low populations.
The entire setup
procedure has now been described. For completeness, here it is in its entirety, together with the declarations at the top of the program. The full NetLogo program is at Modeling Commons.
extensions [ gis table ]
breed [ households household ]
globals [ start-year fertility-age fertility-end-age death-age food-requirement
food-estimate-margin corn-gift-ratio corn-stock-years good-farm-bias
year pdsi-table population-sequence historical-sequence harvest-list fit-badness desired-fit-length ]
patches-own [ value zone patch-quality base-yield historical-base-yield being-farmed? ]
households-own [ harvest age corn-storage unsatisfied-hunger ]
to setup
clear-all ;; as previously described
set-default-shape households "person"
gis-map-load "adata/Map.asc"
set pdsi-table load-zone-file "adata/ZoneAdjPDSI.txt"
ask patches [ initialise-patch ]
set start-year 800 ;; initial year
set food-requirement (5 * 160) ;; a person needs 160 kg of corn each year, while a typical household consist of 5 persons
set fertility-age 17 ;; the minimum age of agents (households) that can reproduce (spin off daughter households)
set fertility-end-age 36 ;; the maximum age of agents (households) that can reproduce
set death-age 38 ;; the age after which agents (households) die (dissolve when the matriarch dies)
set corn-gift-ratio 0.33 ;; each new household gets this fraction of the corn storage of the parent household
set corn-stock-years 2 ;; corn can be stored for this many years
set year start-year
set population-sequence [ ]
set harvest-list [ ]
set historical-sequence
[ 14 14 14 14 14 14 14 14 14 14 13 13 13 12 12 12 12 12 12 12 11 11 11 11 11 10 10 10 10 10 9 9 9 9 9 9 9 9 8 8 7 7 7 7 7 7 7 7 7 7
28 29 29 29 29 29 29 29 29 29 29 30 30 31 31 31 31 32 32 32 32 33 33 33 33 33 37 37 37 37 37 39 39 39 40 40 40 40 41 41 41 42 42
42 42 42 42 42 42 42 56 58 58 58 58 58 58 58 58 58 58 60 60 60 60 60 60 61 61 61 60 61 61 61 61 59 61 61 61 61 61 62 62 62 62 62
62 62 62 62 61 63 63 63 63 63 63 63 63 63 66 67 67 67 67 67 67 67 67 67 66 67 67 67 67 67 67 66 66 66 66 69 69 69 69 65 68 68 68
68 66 67 67 67 66 66 66 66 66 66 66 68 68 68 68 68 68 68 68 68 87 87 87 87 87 87 87 86 86 87 85 86 86 87 87 87 87 88 88 88 87 87
87 87 87 88 90 90 90 90 89 91 92 92 95 95 95 95 97 97 94 94 94 94 95 95 95 95 95 95 134 139 139 139 139 139 139 142 142 142 142
143 143 146 146 146 146 151 151 153 151 151 151 151 151 156 164 164 164 164 163 163 163 163 165 165 165 165 164 164 163 164 164
164 161 162 162 162 162 162 161 164 164 165 165 166 166 166 167 166 166 166 166 159 159 160 160 158 159 160 160 161 162 162 162
148 150 151 151 151 151 151 149 149 147 148 148 148 148 150 150 150 151 151 152 152 152 153 153 153 116 117 118 118 119 119 120
121 122 126 124 124 124 126 127 127 129 131 131 133 132 134 134 136 137 133 138 138 139 139 138 140 140 142 142 142 143 143 144
145 145 147 146 147 147 149 149 150 151 151 145 146 146 147 148 149 149 151 153 154 155 157 159 160 161 163 164 166 167 167 169
169 171 171 173 170 173 176 176 178 178 180 182 184 184 185 188 189 190 192 191 193 192 194 194 194 195 196 199 200 192 193 194
196 196 197 200 202 202 204 201 209 208 211 212 212 213 214 215 216 210 208 206 201 196 188 181 176 172 167 159 156 148 146 141
120 118 114 106 103 95 90 88 83 74 70 68 60 58 56 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 ]
make-region-label 45 98 "Non-arable uplands" black
make-region-label 28 47 "General valley" brown
make-region-label 23 73 "North valley" red
make-region-label 50 68 "Mid valley" gray
make-region-label 25 84 "Dunes" green
make-region-label 20 102 "Kinbiko Canyon" pink
make-region-label 39 11 "Arable uplands" orange
set desired-fit-length (length filter [ ? > 0 ] historical-sequence)
calculate-patch-yields
create-households 14 [
initialise-household (random 29)
set corn-storage (n-values (corn-stock-years + 1) [ 600 + random-float 400 ])
find-best-farm potential-farms
]
reset-ticks
tick-advance start-year
end
The interface for the model looks like this:
Most of the parameters for this model have been given plausible values determined by the literature. The fertility
and usable-farm-fraction
sliders, however, are adjusted to maximise the quality of fit between the model and archaeological reality, and we need to take a closer look at these. The diagram below shows root-mean-square errors for various parameter combinations, averaged over ten runs each. A good value for usable-farm-fraction
is 0.21. A fertility
value of 0.9 generally works well, since it gives a realistically slow population rise, but occasionally it can lead to the population dying out. A value of 0.1 performs almost as well, and is somewhat safer. It is also not far off the value of 0.125 suggested by Axtell et al.
The graph below shows 100 runs of the model (red), compared to the historical data (blue). It can be seen that, in one of the runs, the population dies out early. The other 99 runs follow the historical record reasonably well, however.
This model begs for extension, of course. Not just by restoring the water availability and dwelling location factors I have excluded – the quality of human decision-making could also be improved. The difficulty is that land-use models generally rely on surveying and interviewing farmers in order to reveal their decision processes. The Anasazi are no longer available for such investigation. However, even the simplified model we have presented here has succeeded in replicating the population variation in Long House valley over time. The exception, as Axtell et al. point out, is the final 50-year period. Archaeology shows that the valley was deserted around 1300, but in fact a reasonable population could have survived. Some other factors must have been in operation to cause the mass exodus that occurred. Plausible guesses about those factors could also be included in the model.