Course Plan

This short course introduces a statistical workflow for an experiment here on the NM-AIST campus about the effect of drought, fire, and herbivory on acacia trees. The goal for the week is to learn some shared statistical language so you can work with statisticians (like me) on your research. This week we take the first 4 steps, which help to design an experiment:

  1. What are our scientific goals ?
  2. How will we get data that informs our scientific goals ? What is our experimental design ?
  3. How do we model (i.e. tell a mathematical story about) our data ?
  4. Assuming this model, will we get precise enough estimates to achieve our scientific goals ?

In a future class, after collecting data, we will take the next step:

  1. Is our model realistic enough to achieve our scientific goal ?

Team activity (2 minutes): Throughout class, we will have activities in green. You will work in pairs, and then share your work with the class. A statistical workflow is best done in teams, not alone. Please introduce yourself to your teammate !

Reading

Students are welcome to read these course materials ahead of class. Maoni karibu !

Other reading can include the online textbook Regression and Other Stories, by Gelman et al..

Other references The course materials are based mainly on Regression and Other Stories, by Gelman et al. Chapter 16 about design, Gelman (2023), and Towards A Principled Bayesian Workflow, by Betancourt. For more about step 5 and beyond, see also Bayesian Workflow, by Gelman et al.. We include references below.

Scientific Goals

We want to estimate the effects of drought, fire, and herbivory on African acacia species. Many outcomes are interesting, but in this workflow we focus on height (in centimeters, in September 2023).

Team activity (5 minutes): Why is it important to estimate these effects ? What would you want to know ?

An answer

The height of a tree is a simple way to measure the size and growth rate of a tree. As a tree grows, its height usually increases. Trees that are in poor health or struggling to get key resources will grow very little and thus change little in height. Height may even decrease if upper branches die! So tree growth rate is a useful indicator of the health of a tree.

Tree height is not just a quick way to measure size; the height of a tree itself matters for many reasons. Taller trees can capture more light, which allows them to compete with other plants. Once a tree is tall enough, its leaves are above the height of fire and animals, and so it will escape these stresses. The height of a tree also matters aesthetically; tall trees are beautiful!

Perhaps you can add another reason why tree height is worth measuring?

Acacia Trees

Acacia trees (Miti Migunga) are a very important group of trees across the world, but particularly in Africa. These include such iconic species as the umbrella thorn and the fever tree.

Disney’s Lion King (umbrella thorn acacia ?)
Disney’s Lion King (umbrella thorn acacia ?)


Currently, the former genus Acacia has been divided into two genera: Vachellia and Senegalia.

The twelve Acacia trees in this study are as follows:

Vachellia nilotica, V. gerardii, V. robusta, V. seyal, V. sieberiana, V. tortilis, V. xanthophloea

Senegalia brevispica, S. mellifera, S. polyacantha, S. senegal, V. drepanolobium

These will be referred to using a 6-letter code, created from the first three letters of the genus and species. For instance, Vachellia nilotica is “vacnil” and Senegalia brevispica is “senbre”.

Experimental design

We will get data from a factorial experiment (factors: drought, fire, herbivory) on the NM-AIST campus done by Arjun Potter and Charles Luchagula. They transplanted seedlings of 12 Acacia species from Serengeti National Park to the experimental site on campus.

Arjun Potter collecting a seed pod from Vachellia sieberiana on the NM-AIST campus.
Arjun Potter collecting a seed pod from Vachellia sieberiana on the NM-AIST campus.
Charles Luchagula using a porometer to measure stomatal conductance.
Charles Luchagula using a porometer to measure stomatal conductance.

Layout of the site

The experimental site consists of 10 blocks (12 m × 12 m), 2 m apart. Drought will be implemented via rainout shelters over each block.

Each block contains 9 subplots.

Each subplot contains 6 plant positions.

Team activity (7 minutes): Try to draw this layout on your paper. What questions do you have ?

An answer

Photo of 1 block containing 9 subplots
Photo of 1 block containing 9 subplots

Random assignment of treatment

Fire will be randomly assigned at block-level (5 out of 10 blocks get fire treatment), while drought and herbivory will be randomly assigned at subplot-level.

Team activity (12 minutes): Why do you think the researchers chose to assign fire at the block level and not the subplot level ? What are the advantages and disadvantages ?

An answer

Later in the workflow, we will see there is more precision with random assignment at subplot-level as opposed to block-level. When we assign fire at the block level, we cannot distinguish between the effect of being in a block versus the effect of fire. For example, if block 3 got fire but block 7 did not get fire, is the difference between the plant heights because of the fire treatment ? Or because block 3 had slightly different soil than block 7 ? Or maybe their rainout shelters are built slightly differently ? When treatment (herbivory and drought) is assigned at a lower level like a subplot, we can compare within the same block and isolate the effect of treatment alone.

But manipulating fire is more difficult than manipulating drought or herbivory. It would be difficult to do the controlled burns at the subplot level.


Of the 9 subplots per block, one is reserved for a corn seedling to capture environment differences across blocks (corns are clones of each other, unlike acacias where there is much genetic variation even within species). The remaining 8 subplots are randomly assigned to 4 treatment combinations: control, drought, herbivory, drought + herbivory. Each pair of subplots (6 + 6 plant positions) with the same treatment combination gets all 12 Acacia species, arranged randomly among plant positions.

Team activity (8 minutes): Focus on one block, with its 9 subplots. Draw one possible randomization on your paper, which subplots get which treatments ? Where are the 12 Acacia species planted (you can label then 1 through 12) ?

An answer

Implementation of treatment

Fire will be implemented via dried grasses lit with a drip torch on four sides of the fire-assigned blocks.

Drought will be implemented via rainout shelters over each block. All rainfall will be diverted into ditches. Garden hoses will be used to water the subplots. Control subplots will be watered regularly on a set schedule; the amount of water added will be comparable to average rainfall conditions in East African savannas. Drought subplots will receive no water for a period of weeks to months. This mimics the rainfall uncertainty that is common in East Africa. This is an el Niño year. Dry and wet seasons have been getting more chaotic. For this reason, the experiment does not try to match the rainfall at Nelson Mandela. There are 2 soil moisture probes in each rainout shelter. These will be used to quantify exactly how dry the soil gets. Even when soil looks or feels dry, there is actually a small amount of water bound to the soil. Plants that grow in dry regions can suck this tiny amount of water from the soil !

Herbivory will be implemented via clipping and stripping leaves.

Team activity (9 minutes): Do you have any implementation questions or concerns ?

An answer

A control (watered) subplot might be next to a drought (no water) subplot. If the drought subplot could access the water from the control subplot, its treatment becomes more complicated. This interference between experimental units is discussed in “Causal Inference: What If” book by Hernán and Robins Fine Point 1.1 on page 5.

Because the soil is a heavy clay, water does not flow quickly through the soil. For this reason, small berms and trenches can stop water from flowing between the control (watered) and drought (no water) subplots.

Mathematical notation

To clarify our design and data, we define mathematical notation.

  • \(i =\) a tree, a row in our data.

Experimental design:

  • treatment factors:
    • \(\text{Fire}_i = 1\) if tree \(i\) gets fire treatment, \(=0\) otherwise
    • \(\text{Drought}_i = 1\) if tree \(i\) gets drought treatment, \(=0\) otherwise
    • \(\text{Herbivory}_i = 1\) if tree \(i\) gets herbivory treatment, \(=0\) otherwise
  • \(s[i] =\) species of tree \(i\), one of 12 species of Acacia being studied. The names of the trees are abbreviated; for instance, “vacxan” is Vachellia xanthophloea, known as “Mgunga Maji” in Kiswahili. Each species is given a unique number (Species_number).
  • tree location:
    • \(b[i] =\) block of tree \(i\), one of 10 blocks
    • \(sp[i] =\) subplot of tree \(i\), one of 9 subplots nested in a block
    • \(p[i] =\) plant position of tree \(i\), one of 6 positions nested in a subplot
  • Data:
    • \(y_i =\) height (in centimeters, in October 2023) of tree \(i\), not yet observed
    • \(\text{corn}_b =\) height of corn in block \(b\) (in centimeters, in October 2023), not yet observed

Here is the experimental design. Each row is a tree, columns are variables: location, species, and treatment factors).

tree Fire Drought Herbivory Species Species_number Block Subplot PlantPosition is_corn
1 0 1 1 vacsey 9 1 1 1 0
2 0 1 1 senpol 4 1 1 2 0
3 0 1 1 vacger 6 1 1 3 0
4 0 1 1 albhar 1 1 1 4 0
5 0 1 1 senbre 2 1 1 5 0
6 0 1 1 sensen 5 1 1 6 0
7 0 0 1 albhar 1 1 2 1 0
8 0 0 1 vacrob 8 1 2 2 0
9 0 0 1 senmel 3 1 2 3 0
10 0 0 1 vactor 11 1 2 4 0
11 0 0 1 vacxan 12 1 2 5 0
12 0 0 1 vacnil 7 1 2 6 0
13 0 0 0 albhar 1 1 3 1 0
14 0 0 0 vacsie 10 1 3 2 0
15 0 0 0 vacrob 8 1 3 3 0
16 0 0 0 senmel 3 1 3 4 0
17 0 0 0 senbre 2 1 3 5 0
18 0 0 0 sensen 5 1 3 6 0
19 0 1 0 vacxan 12 1 4 1 0
20 0 1 0 vactor 11 1 4 2 0
21 0 1 0 vacrob 8 1 4 3 0
22 0 1 0 senpol 4 1 4 4 0
23 0 1 0 senbre 2 1 4 5 0
24 0 1 0 senmel 3 1 4 6 0
25 0 0 1 sensen 5 1 5 1 0
26 0 0 1 senbre 2 1 5 2 0
27 0 0 1 vacsie 10 1 5 3 0
28 0 0 1 vacsey 9 1 5 4 0
29 0 0 1 senpol 4 1 5 5 0
30 0 0 1 vacger 6 1 5 6 0
31 0 1 1 senmel 3 1 6 1 0
32 0 1 1 vacrob 8 1 6 2 0
33 0 1 1 vacnil 7 1 6 3 0
34 0 1 1 vactor 11 1 6 4 0
35 0 1 1 vacxan 12 1 6 5 0
36 0 1 1 vacsie 10 1 6 6 0
37 0 0 0 vacger 6 1 7 1 1
38 0 0 0 vacnil 7 1 8 1 0
39 0 0 0 senpol 4 1 8 2 0
40 0 0 0 vactor 11 1 8 3 0
41 0 0 0 vacsey 9 1 8 4 0
42 0 0 0 vacger 6 1 8 5 0
43 0 0 0 vacxan 12 1 8 6 0
44 0 1 0 sensen 5 1 9 1 0
45 0 1 0 vacsie 10 1 9 2 0
46 0 1 0 vacnil 7 1 9 3 0
47 0 1 0 vacger 6 1 9 4 0
48 0 1 0 albhar 1 1 9 5 0
49 0 1 0 vacsey 9 1 9 6 0
50 1 0 0 vacrob 8 2 1 1 1
51 1 0 0 vacxan 12 2 2 1 0
52 1 0 0 senbre 2 2 2 2 0
53 1 0 0 sensen 5 2 2 3 0
54 1 0 0 vacrob 8 2 2 4 0
55 1 0 0 senpol 4 2 2 5 0
56 1 0 0 vactor 11 2 2 6 0
57 1 0 1 vacsey 9 2 3 1 0
58 1 0 1 vacger 6 2 3 2 0
59 1 0 1 senpol 4 2 3 3 0
60 1 0 1 vactor 11 2 3 4 0
61 1 0 1 sensen 5 2 3 5 0
62 1 0 1 vacrob 8 2 3 6 0
63 1 1 1 vacsie 10 2 4 1 0
64 1 1 1 vacger 6 2 4 2 0
65 1 1 1 senmel 3 2 4 3 0
66 1 1 1 vacnil 7 2 4 4 0
67 1 1 1 vacrob 8 2 4 5 0
68 1 1 1 vacsey 9 2 4 6 0
69 1 1 0 albhar 1 2 5 1 0
70 1 1 0 senmel 3 2 5 2 0
71 1 1 0 vactor 11 2 5 3 0
72 1 1 0 vacrob 8 2 5 4 0
73 1 1 0 vacnil 7 2 5 5 0
74 1 1 0 vacger 6 2 5 6 0
75 1 1 0 senpol 4 2 6 1 0
76 1 1 0 vacsie 10 2 6 2 0
77 1 1 0 vacxan 12 2 6 3 0
78 1 1 0 vacsey 9 2 6 4 0
79 1 1 0 senbre 2 2 6 5 0
80 1 1 0 sensen 5 2 6 6 0
81 1 0 0 vacnil 7 2 7 1 0
82 1 0 0 albhar 1 2 7 2 0
83 1 0 0 vacsey 9 2 7 3 0
84 1 0 0 vacger 6 2 7 4 0
85 1 0 0 vacsie 10 2 7 5 0
86 1 0 0 senmel 3 2 7 6 0
87 1 0 1 vacsie 10 2 8 1 0
88 1 0 1 vacxan 12 2 8 2 0
89 1 0 1 vacnil 7 2 8 3 0
90 1 0 1 albhar 1 2 8 4 0
91 1 0 1 senmel 3 2 8 5 0
92 1 0 1 senbre 2 2 8 6 0
93 1 1 1 senbre 2 2 9 1 0
94 1 1 1 sensen 5 2 9 2 0
95 1 1 1 vactor 11 2 9 3 0
96 1 1 1 senpol 4 2 9 4 0
97 1 1 1 vacxan 12 2 9 5 0
98 1 1 1 albhar 1 2 9 6 0
99 1 0 0 vacnil 7 3 1 1 0
100 1 0 0 vacsie 10 3 1 2 0
101 1 0 0 vacger 6 3 1 3 0
102 1 0 0 vactor 11 3 1 4 0
103 1 0 0 sensen 5 3 1 5 0
104 1 0 0 vacrob 8 3 1 6 0
105 1 0 1 albhar 1 3 2 1 0
106 1 0 1 vacnil 7 3 2 2 0
107 1 0 1 vactor 11 3 2 3 0
108 1 0 1 senpol 4 3 2 4 0
109 1 0 1 vacxan 12 3 2 5 0
110 1 0 1 sensen 5 3 2 6 0
111 1 1 1 vacsie 10 3 3 1 0
112 1 1 1 vacger 6 3 3 2 0
113 1 1 1 senmel 3 3 3 3 0
114 1 1 1 senbre 2 3 3 4 0
115 1 1 1 vacsey 9 3 3 5 0
116 1 1 1 albhar 1 3 3 6 0
117 1 0 0 senbre 2 3 4 1 1
118 1 1 0 senpol 4 3 5 1 0
119 1 1 0 albhar 1 3 5 2 0
120 1 1 0 senmel 3 3 5 3 0
121 1 1 0 senbre 2 3 5 4 0
122 1 1 0 vacsey 9 3 5 5 0
123 1 1 0 vacnil 7 3 5 6 0
124 1 1 1 vacxan 12 3 6 1 0
125 1 1 1 vacnil 7 3 6 2 0
126 1 1 1 vacrob 8 3 6 3 0
127 1 1 1 sensen 5 3 6 4 0
128 1 1 1 senpol 4 3 6 5 0
129 1 1 1 vactor 11 3 6 6 0
130 1 1 0 vacrob 8 3 7 1 0
131 1 1 0 vacsie 10 3 7 2 0
132 1 1 0 sensen 5 3 7 3 0
133 1 1 0 vacger 6 3 7 4 0
134 1 1 0 vacxan 12 3 7 5 0
135 1 1 0 vactor 11 3 7 6 0
136 1 0 1 vacger 6 3 8 1 0
137 1 0 1 vacrob 8 3 8 2 0
138 1 0 1 senmel 3 3 8 3 0
139 1 0 1 senbre 2 3 8 4 0
140 1 0 1 vacsie 10 3 8 5 0
141 1 0 1 vacsey 9 3 8 6 0
142 1 0 0 senmel 3 3 9 1 0
143 1 0 0 vacxan 12 3 9 2 0
144 1 0 0 albhar 1 3 9 3 0
145 1 0 0 vacsey 9 3 9 4 0
146 1 0 0 senpol 4 3 9 5 0
147 1 0 0 senbre 2 3 9 6 0
148 0 0 0 vacxan 12 4 1 1 1
149 0 1 0 vacsey 9 4 2 1 0
150 0 1 0 vacsie 10 4 2 2 0
151 0 1 0 vactor 11 4 2 3 0
152 0 1 0 senbre 2 4 2 4 0
153 0 1 0 vacxan 12 4 2 5 0
154 0 1 0 senmel 3 4 2 6 0
155 0 0 0 sensen 5 4 3 1 0
156 0 0 0 vacger 6 4 3 2 0
157 0 0 0 vacxan 12 4 3 3 0
158 0 0 0 vacnil 7 4 3 4 0
159 0 0 0 senpol 4 4 3 5 0
160 0 0 0 vactor 11 4 3 6 0
161 0 1 1 vacger 6 4 4 1 0
162 0 1 1 senpol 4 4 4 2 0
163 0 1 1 vacsey 9 4 4 3 0
164 0 1 1 vacrob 8 4 4 4 0
165 0 1 1 senmel 3 4 4 5 0
166 0 1 1 sensen 5 4 4 6 0
167 0 1 0 albhar 1 4 5 1 0
168 0 1 0 vacrob 8 4 5 2 0
169 0 1 0 vacger 6 4 5 3 0
170 0 1 0 senpol 4 4 5 4 0
171 0 1 0 vacnil 7 4 5 5 0
172 0 1 0 sensen 5 4 5 6 0
173 0 0 1 senmel 3 4 6 1 0
174 0 0 1 vactor 11 4 6 2 0
175 0 0 1 vacrob 8 4 6 3 0
176 0 0 1 vacsey 9 4 6 4 0
177 0 0 1 senbre 2 4 6 5 0
178 0 0 1 senpol 4 4 6 6 0
179 0 0 1 vacsie 10 4 7 1 0
180 0 0 1 vacnil 7 4 7 2 0
181 0 0 1 vacxan 12 4 7 3 0
182 0 0 1 albhar 1 4 7 4 0
183 0 0 1 vacger 6 4 7 5 0
184 0 0 1 sensen 5 4 7 6 0
185 0 0 0 senbre 2 4 8 1 0
186 0 0 0 albhar 1 4 8 2 0
187 0 0 0 vacrob 8 4 8 3 0
188 0 0 0 senmel 3 4 8 4 0
189 0 0 0 vacsey 9 4 8 5 0
190 0 0 0 vacsie 10 4 8 6 0
191 0 1 1 vacsie 10 4 9 1 0
192 0 1 1 vactor 11 4 9 2 0
193 0 1 1 senbre 2 4 9 3 0
194 0 1 1 albhar 1 4 9 4 0
195 0 1 1 vacxan 12 4 9 5 0
196 0 1 1 vacnil 7 4 9 6 0
197 0 0 0 vacsey 9 5 1 1 0
198 0 0 0 vacxan 12 5 1 2 0
199 0 0 0 sensen 5 5 1 3 0
200 0 0 0 senmel 3 5 1 4 0
201 0 0 0 vacnil 7 5 1 5 0
202 0 0 0 vactor 11 5 1 6 0
203 0 0 1 vacnil 7 5 2 1 0
204 0 0 1 senbre 2 5 2 2 0
205 0 0 1 vactor 11 5 2 3 0
206 0 0 1 vacrob 8 5 2 4 0
207 0 0 1 vacxan 12 5 2 5 0
208 0 0 1 vacsey 9 5 2 6 0
209 0 1 0 vacrob 8 5 3 1 0
210 0 1 0 vactor 11 5 3 2 0
211 0 1 0 vacsey 9 5 3 3 0
212 0 1 0 vacsie 10 5 3 4 0
213 0 1 0 senmel 3 5 3 5 0
214 0 1 0 vacger 6 5 3 6 0
215 0 0 0 vacrob 8 5 4 1 0
216 0 0 0 vacger 6 5 4 2 0
217 0 0 0 senpol 4 5 4 3 0
218 0 0 0 vacsie 10 5 4 4 0
219 0 0 0 albhar 1 5 4 5 0
220 0 0 0 senbre 2 5 4 6 0
221 0 1 1 vacnil 7 5 5 1 0
222 0 1 1 albhar 1 5 5 2 0
223 0 1 1 senpol 4 5 5 3 0
224 0 1 1 sensen 5 5 5 4 0
225 0 1 1 vactor 11 5 5 5 0
226 0 1 1 senmel 3 5 5 6 0
227 0 0 1 senpol 4 5 6 1 0
228 0 0 1 vacsie 10 5 6 2 0
229 0 0 1 vacger 6 5 6 3 0
230 0 0 1 albhar 1 5 6 4 0
231 0 0 1 senmel 3 5 6 5 0
232 0 0 1 sensen 5 5 6 6 0
233 0 1 1 vacsie 10 5 7 1 0
234 0 1 1 vacsey 9 5 7 2 0
235 0 1 1 vacxan 12 5 7 3 0
236 0 1 1 vacrob 8 5 7 4 0
237 0 1 1 vacger 6 5 7 5 0
238 0 1 1 senbre 2 5 7 6 0
239 0 0 0 vactor 11 5 8 1 1
240 0 1 0 sensen 5 5 9 1 0
241 0 1 0 albhar 1 5 9 2 0
242 0 1 0 vacxan 12 5 9 3 0
243 0 1 0 senpol 4 5 9 4 0
244 0 1 0 senbre 2 5 9 5 0
245 0 1 0 vacnil 7 5 9 6 0
246 0 0 0 sensen 5 6 1 1 0
247 0 0 0 albhar 1 6 1 2 0
248 0 0 0 senbre 2 6 1 3 0
249 0 0 0 vactor 11 6 1 4 0
250 0 0 0 vacrob 8 6 1 5 0
251 0 0 0 vacxan 12 6 1 6 0
252 0 0 1 vacxan 12 6 2 1 0
253 0 0 1 vacsie 10 6 2 2 0
254 0 0 1 senbre 2 6 2 3 0
255 0 0 1 albhar 1 6 2 4 0
256 0 0 1 vacsey 9 6 2 5 0
257 0 0 1 senpol 4 6 2 6 0
258 0 1 1 vacsey 9 6 3 1 0
259 0 1 1 vacnil 7 6 3 2 0
260 0 1 1 vactor 11 6 3 3 0
261 0 1 1 senpol 4 6 3 4 0
262 0 1 1 senmel 3 6 3 5 0
263 0 1 1 sensen 5 6 3 6 0
264 0 1 0 sensen 5 6 4 1 0
265 0 1 0 vacger 6 6 4 2 0
266 0 1 0 senpol 4 6 4 3 0
267 0 1 0 vacsie 10 6 4 4 0
268 0 1 0 vacxan 12 6 4 5 0
269 0 1 0 vactor 11 6 4 6 0
270 0 0 0 senmel 3 6 5 1 1
271 0 0 1 vacger 6 6 6 1 0
272 0 0 1 sensen 5 6 6 2 0
273 0 0 1 vacrob 8 6 6 3 0
274 0 0 1 vactor 11 6 6 4 0
275 0 0 1 vacnil 7 6 6 5 0
276 0 0 1 senmel 3 6 6 6 0
277 0 1 0 senbre 2 6 7 1 0
278 0 1 0 vacsey 9 6 7 2 0
279 0 1 0 vacrob 8 6 7 3 0
280 0 1 0 albhar 1 6 7 4 0
281 0 1 0 vacnil 7 6 7 5 0
282 0 1 0 senmel 3 6 7 6 0
283 0 1 1 vacrob 8 6 8 1 0
284 0 1 1 vacger 6 6 8 2 0
285 0 1 1 albhar 1 6 8 3 0
286 0 1 1 vacxan 12 6 8 4 0
287 0 1 1 vacsie 10 6 8 5 0
288 0 1 1 senbre 2 6 8 6 0
289 0 0 0 vacsie 10 6 9 1 0
290 0 0 0 senmel 3 6 9 2 0
291 0 0 0 vacger 6 6 9 3 0
292 0 0 0 vacnil 7 6 9 4 0
293 0 0 0 senpol 4 6 9 5 0
294 0 0 0 vacsey 9 6 9 6 0
295 1 1 0 vacsey 9 7 1 1 0
296 1 1 0 senmel 3 7 1 2 0
297 1 1 0 vacrob 8 7 1 3 0
298 1 1 0 vacger 6 7 1 4 0
299 1 1 0 senpol 4 7 1 5 0
300 1 1 0 vacnil 7 7 1 6 0
301 1 0 0 senmel 3 7 2 1 0
302 1 0 0 vacsie 10 7 2 2 0
303 1 0 0 vacnil 7 7 2 3 0
304 1 0 0 vacger 6 7 2 4 0
305 1 0 0 albhar 1 7 2 5 0
306 1 0 0 vacrob 8 7 2 6 0
307 1 0 1 albhar 1 7 3 1 0
308 1 0 1 vactor 11 7 3 2 0
309 1 0 1 vacsey 9 7 3 3 0
310 1 0 1 vacger 6 7 3 4 0
311 1 0 1 vacsie 10 7 3 5 0
312 1 0 1 vacnil 7 7 3 6 0
313 1 0 0 vactor 11 7 4 1 0
314 1 0 0 vacsey 9 7 4 2 0
315 1 0 0 vacxan 12 7 4 3 0
316 1 0 0 senpol 4 7 4 4 0
317 1 0 0 sensen 5 7 4 5 0
318 1 0 0 senbre 2 7 4 6 0
319 1 1 1 senpol 4 7 5 1 0
320 1 1 1 senbre 2 7 5 2 0
321 1 1 1 sensen 5 7 5 3 0
322 1 1 1 vactor 11 7 5 4 0
323 1 1 1 vacxan 12 7 5 5 0
324 1 1 1 senmel 3 7 5 6 0
325 1 1 1 vacger 6 7 6 1 0
326 1 1 1 albhar 1 7 6 2 0
327 1 1 1 vacnil 7 7 6 3 0
328 1 1 1 vacrob 8 7 6 4 0
329 1 1 1 vacsie 10 7 6 5 0
330 1 1 1 vacsey 9 7 6 6 0
331 1 0 1 sensen 5 7 7 1 0
332 1 0 1 senpol 4 7 7 2 0
333 1 0 1 vacrob 8 7 7 3 0
334 1 0 1 senbre 2 7 7 4 0
335 1 0 1 vacxan 12 7 7 5 0
336 1 0 1 senmel 3 7 7 6 0
337 1 0 0 albhar 1 7 8 1 1
338 1 1 0 sensen 5 7 9 1 0
339 1 1 0 vactor 11 7 9 2 0
340 1 1 0 albhar 1 7 9 3 0
341 1 1 0 vacsie 10 7 9 4 0
342 1 1 0 senbre 2 7 9 5 0
343 1 1 0 vacxan 12 7 9 6 0
344 1 1 0 vacger 6 8 1 1 0
345 1 1 0 albhar 1 8 1 2 0
346 1 1 0 senpol 4 8 1 3 0
347 1 1 0 vacsie 10 8 1 4 0
348 1 1 0 senbre 2 8 1 5 0
349 1 1 0 sensen 5 8 1 6 0
350 1 0 0 sensen 5 8 2 1 0
351 1 0 0 senmel 3 8 2 2 0
352 1 0 0 senpol 4 8 2 3 0
353 1 0 0 vacsey 9 8 2 4 0
354 1 0 0 vacger 6 8 2 5 0
355 1 0 0 vactor 11 8 2 6 0
356 1 0 0 vacsey 9 8 3 1 1
357 1 1 0 vacsey 9 8 4 1 0
358 1 1 0 vactor 11 8 4 2 0
359 1 1 0 vacnil 7 8 4 3 0
360 1 1 0 vacrob 8 8 4 4 0
361 1 1 0 senmel 3 8 4 5 0
362 1 1 0 vacxan 12 8 4 6 0
363 1 0 1 senbre 2 8 5 1 0
364 1 0 1 senpol 4 8 5 2 0
365 1 0 1 senmel 3 8 5 3 0
366 1 0 1 vactor 11 8 5 4 0
367 1 0 1 vacger 6 8 5 5 0
368 1 0 1 vacsie 10 8 5 6 0
369 1 1 1 senbre 2 8 6 1 0
370 1 1 1 sensen 5 8 6 2 0
371 1 1 1 vacsie 10 8 6 3 0
372 1 1 1 vacrob 8 8 6 4 0
373 1 1 1 vacxan 12 8 6 5 0
374 1 1 1 senpol 4 8 6 6 0
375 1 0 1 vacxan 12 8 7 1 0
376 1 0 1 albhar 1 8 7 2 0
377 1 0 1 vacrob 8 8 7 3 0
378 1 0 1 sensen 5 8 7 4 0
379 1 0 1 vacnil 7 8 7 5 0
380 1 0 1 vacsey 9 8 7 6 0
381 1 1 1 albhar 1 8 8 1 0
382 1 1 1 vacnil 7 8 8 2 0
383 1 1 1 vacger 6 8 8 3 0
384 1 1 1 vactor 11 8 8 4 0
385 1 1 1 senmel 3 8 8 5 0
386 1 1 1 vacsey 9 8 8 6 0
387 1 0 0 vacrob 8 8 9 1 0
388 1 0 0 vacsie 10 8 9 2 0
389 1 0 0 senbre 2 8 9 3 0
390 1 0 0 albhar 1 8 9 4 0
391 1 0 0 vacnil 7 8 9 5 0
392 1 0 0 vacxan 12 8 9 6 0
393 0 1 0 senmel 3 9 1 1 0
394 0 1 0 vacxan 12 9 1 2 0
395 0 1 0 albhar 1 9 1 3 0
396 0 1 0 senbre 2 9 1 4 0
397 0 1 0 vacrob 8 9 1 5 0
398 0 1 0 vacnil 7 9 1 6 0
399 0 0 1 vacsey 9 9 2 1 0
400 0 0 1 sensen 5 9 2 2 0
401 0 0 1 albhar 1 9 2 3 0
402 0 0 1 senpol 4 9 2 4 0
403 0 0 1 senmel 3 9 2 5 0
404 0 0 1 vacger 6 9 2 6 0
405 0 0 0 sensen 5 9 3 1 1
406 0 0 0 sensen 5 9 4 1 0
407 0 0 0 vacsie 10 9 4 2 0
408 0 0 0 vacrob 8 9 4 3 0
409 0 0 0 senbre 2 9 4 4 0
410 0 0 0 albhar 1 9 4 5 0
411 0 0 0 vactor 11 9 4 6 0
412 0 1 1 vacxan 12 9 5 1 0
413 0 1 1 albhar 1 9 5 2 0
414 0 1 1 sensen 5 9 5 3 0
415 0 1 1 vacnil 7 9 5 4 0
416 0 1 1 vacsey 9 9 5 5 0
417 0 1 1 vactor 11 9 5 6 0
418 0 1 1 senmel 3 9 6 1 0
419 0 1 1 senbre 2 9 6 2 0
420 0 1 1 vacger 6 9 6 3 0
421 0 1 1 senpol 4 9 6 4 0
422 0 1 1 vacrob 8 9 6 5 0
423 0 1 1 vacsie 10 9 6 6 0
424 0 0 0 vacger 6 9 7 1 0
425 0 0 0 senpol 4 9 7 2 0
426 0 0 0 vacsey 9 9 7 3 0
427 0 0 0 vacxan 12 9 7 4 0
428 0 0 0 senmel 3 9 7 5 0
429 0 0 0 vacnil 7 9 7 6 0
430 0 1 0 senpol 4 9 8 1 0
431 0 1 0 vacsey 9 9 8 2 0
432 0 1 0 vactor 11 9 8 3 0
433 0 1 0 sensen 5 9 8 4 0
434 0 1 0 vacsie 10 9 8 5 0
435 0 1 0 vacger 6 9 8 6 0
436 0 0 1 vactor 11 9 9 1 0
437 0 0 1 vacxan 12 9 9 2 0
438 0 0 1 vacsie 10 9 9 3 0
439 0 0 1 senbre 2 9 9 4 0
440 0 0 1 vacrob 8 9 9 5 0
441 0 0 1 vacnil 7 9 9 6 0
442 1 1 0 senmel 3 10 1 1 0
443 1 1 0 vactor 11 10 1 2 0
444 1 1 0 albhar 1 10 1 3 0
445 1 1 0 vacsie 10 10 1 4 0
446 1 1 0 vacxan 12 10 1 5 0
447 1 1 0 vacger 6 10 1 6 0
448 1 0 0 senpol 4 10 2 1 0
449 1 0 0 senmel 3 10 2 2 0
450 1 0 0 sensen 5 10 2 3 0
451 1 0 0 senbre 2 10 2 4 0
452 1 0 0 vacsie 10 10 2 5 0
453 1 0 0 vacrob 8 10 2 6 0
454 1 0 1 vacrob 8 10 3 1 0
455 1 0 1 albhar 1 10 3 2 0
456 1 0 1 senmel 3 10 3 3 0
457 1 0 1 vactor 11 10 3 4 0
458 1 0 1 vacsie 10 10 3 5 0
459 1 0 1 vacnil 7 10 3 6 0
460 1 0 0 vacnil 7 10 4 1 1
461 1 0 1 vacger 6 10 5 1 0
462 1 0 1 vacxan 12 10 5 2 0
463 1 0 1 senpol 4 10 5 3 0
464 1 0 1 sensen 5 10 5 4 0
465 1 0 1 vacsey 9 10 5 5 0
466 1 0 1 senbre 2 10 5 6 0
467 1 1 0 sensen 5 10 6 1 0
468 1 1 0 vacrob 8 10 6 2 0
469 1 1 0 vacsey 9 10 6 3 0
470 1 1 0 senpol 4 10 6 4 0
471 1 1 0 senbre 2 10 6 5 0
472 1 1 0 vacnil 7 10 6 6 0
473 1 1 1 senbre 2 10 7 1 0
474 1 1 1 albhar 1 10 7 2 0
475 1 1 1 vactor 11 10 7 3 0
476 1 1 1 vacsey 9 10 7 4 0
477 1 1 1 vacnil 7 10 7 5 0
478 1 1 1 senmel 3 10 7 6 0
479 1 1 1 vacsie 10 10 8 1 0
480 1 1 1 senpol 4 10 8 2 0
481 1 1 1 vacger 6 10 8 3 0
482 1 1 1 vacrob 8 10 8 4 0
483 1 1 1 vacxan 12 10 8 5 0
484 1 1 1 sensen 5 10 8 6 0
485 1 0 0 vacxan 12 10 9 1 0
486 1 0 0 vactor 11 10 9 2 0
487 1 0 0 vacsey 9 10 9 3 0
488 1 0 0 vacnil 7 10 9 4 0
489 1 0 0 vacger 6 10 9 5 0
490 1 0 0 albhar 1 10 9 6 0

Team activity (6 minutes): what is \(s[3]\) ?

An answer

Tree \(i=3\) in row 3 of the data is the “vacger” species. If you look at the Species_number column, you can see that I’ve labeled this Species_number 6 (out of 12 species we are studying).

So,

\(s[3] = 6\).

What this means is the species number for the third row (the third plant) is number 6.


Modeling

To estimate the effects of drought, fire, and herbivory on the acacia species, we compare heights across groups.

But we have small numbers of trees in each group.

Team activity (8 minutes): How many total trees do we have ? How many trees are Vachellia gerardii that got fire, drought, and not herbivory ?

An answer

10 blocks * (1 corn plant + 8 regular subplots * 6 plants) = 490 total trees.

5 out of the 10 blocks got fire. In each of these 5 blocks, there is a pair of subplots that got drought and not herbivory. But the species Vachellia gerardii is only in one of them, since the 12 acacia species we are studying are split across the pair of subplots. So there are only 5 Vachellia gerardii plants that got fire, drought, and not herbivory.


We could use similar groups (Vachellia gerardii plants that got fire, drought, and herbivory) to help get better estimates of effects.

Team activity (9 minutes): We want to study drought, herbivory, and fire. What are some things we don’t want to study that might affect tree height ?

An answer

We may be able to explain some of the variation in heights by a plant’s location: how far south or east is it ? where is it within a block’s rainout shelter ? And maybe some rainout shelters are slightly different than others in a way that affects plant height.


We combine data across groups and our knowledge of the experimental design together into a statistical model, a mathematical story of how the data could be generated. We use probability distributions to express both randomness in nature and our uncertainty. For example, a block’s rainout shelter could randomly tilt in the wind, which we express with block effects drawn from a probability distribution. As another example, we have uncertainty about the effects of fire, so we express the fire effects as being drawn from a probability distribution.

A useful statistical model will help us get better estimates than simply comparing heights across groups.

The big picture

A mathematical story of how tree height could be generated might include treatment, species, location, and the specific tree:

\[\begin{align*} y_i &= \text{average height} \\ &\qquad + \text{factor (fire, drought, herbivory) effects} \\ &\qquad + \text{species-specific effects} \\ &\qquad + \text{location effects} \\ &\qquad + \text{effect for tree }i \end{align*}\]

Team activity (6 minutes): Which parts of the model do you think will be most important to predict acacia tree height ?

An answer

I don’t have an answer here, because I don’t know about acacia trees. But as a statistician, I imagine that the goal of the experiment is for location effects to be minimized so that we can focus on the effects of interest: fire, drought, herbivory, and variation across species.


The above model is not fully in math yet. Next we will use probability distributions to express each part mathematically. Above, we defined mathematical notation for our data using English letters (like \(y_i\) and \(s[i]\)). Mathematical notation for the model usually uses Greek letters (like \(\beta\) and \(\epsilon\)).

Simulating fake data from the model

As we build our mathematical story (the model) about tree heights, we want to see if it is reasonable. We already have prior information about acacia trees: how tall do we expect them to get after a few months ? what height would be surprising ? how do we expect them to react to fire, drought, and herbivory ? how much do we expect different species to be different ? what differences would be surprising ?

To see if the model agrees with our prior information, we will simulate fake data from the model. This way we can check to see if it feels reasonable to us.

If the fake data seem reasonable enough, we can even analyze that fake data to see how precise we expect our results to be.

Let’s now go one step at a time through our mathematical story.

Average height

Team activity (9 minutes): What might be the average height of the trees across all species for the control (no fire, no drought, no herbivory) ? Draw a distribution to express which you think are the most likely heights. This is your prior about average acacia tree heights.

An answer

Maybe 20 centimeters is most likely ? And probably not less than 10 centimeters or more than 30 centimeters ?

We can express this mathematically with a Normal distribution with mean 20 and standard deviation 5:

\(\beta_0 \sim N(20,5)\)

We will simulate from this distribution to create fake data from our model’s story:

set.seed(1)
beta_0 = rnorm(n = 1, mean = 20, sd = 5)


Factor effects

Let’s write a small part our model for the tree heights in math:

\[\begin{align*} y_i &= \beta_0 + \beta^{\text{fire}} \text{Fire}_i + \beta^{\text{drought}} \text{Drought}_i + \beta^{\text{herbivory}} \text{Herbivory}_i + \epsilon_i \end{align*}\]

Recall our notation:

  • \(i =\) a tree
  • \(\text{Fire}_i = 1\) if tree \(i\) gets fire treatment, \(=0\) otherwise
  • \(\text{Drought}_i = 1\) if tree \(i\) gets drought treatment, \(=0\) otherwise
  • \(\text{Herbivory}_i = 1\) if tree \(i\) gets herbivory treatment, \(=0\) otherwise
  • \(y_i =\) height of tree \(i\)

We also add \(\epsilon_i\), an “error” term that captures all the ways tree \(i\) is special that are not in our model’s prediction \(\beta_0 + \beta^{\text{fire}} \text{Fire}_i + \beta^{\text{drought}} \text{Drought}_i + \beta^{\text{herbivory}} \text{Herbivory}_i\).

Team activity (9 minutes): What does our model predict will be the height of a tree in the control group ?

An answer

The control group gets none of the treatments, so \(\text{Fire}_i = 0\), \(\text{Drought}_i = 0\), and \(\text{Herbivory}_i = 0\). So the model predicts the control height will be:

\[\beta_0 + \beta^{\text{fire}}*0 + \beta^{\text{drought}}*0 + \beta^{\text{herbivory}}*0\] which is just \(\beta_0\).

Team activity (9 minutes): What does our model predict will be the height of a tree that gets fire treatment but not drought or herbivory ?

An answer

This group gets \(\text{Fire}_i = 1\), but \(\text{Drought}_i = 0\) and \(\text{Herbivory}_i = 0\). So the model predicts their height will be:

\[\beta_0 + \beta^{\text{fire}}*1 + \beta^{\text{drought}}*0 + \beta^{\text{herbivory}}*0\] which is \(\beta_0 + \beta^{\text{fire}}\).

Team activity (9 minutes): What does our model predict will be the height of a tree that gets fire and drought but not the herbivory treatment ?

An answer

This group gets \(\text{Fire}_i = 1\) and \(\text{Drought}_i = 1\), but \(\text{Herbivory}_i = 0\). So the model predicts their height will be:

\[\beta_0 + \beta^{\text{fire}}*1 + \beta^{\text{drought}}*1 + \beta^{\text{herbivory}}*0\] which is \(\beta_0 + \beta^{\text{fire}} + \beta^{\text{drought}}\).

Team activity (9 minutes): What might be the average effect of fire (in the absence of drought and herbivory) ? What is your best guess ? Reducing height by 1 centimeter ? 5 centimeters ? Could fire ever be helpful for growth ? Express your guesses mathematically with a Normal distribution. This is your prior about the effect of fire.

An answer

Maybe fire reduces height by 2 centimeters ? Probably it is never helpful and never worse than reducing height by 4 centimeters ?

We can express this mathematically with a Normal distribution with mean -2 and standard deviation 1:

\(\beta^{\text{fire}} \sim N(-2,1)\)

Again, we simulate to create fake data from our model’s story:

beta_fire = rnorm(n = 1, mean = -2, sd = 1)

Team activity (11 minutes): Repeat the above for the effects of drought and herbivory.

An answer

Similarly, let’s think about the effects of drought and herbivory:

\(\beta^{\text{drought}} \sim N(-4,2)\)

beta_drought = rnorm(n = 1, mean = -4, sd = 2)

\(\beta^{\text{herbivory}} \sim N(-2,1)\)

beta_herbivory = rnorm(n = 1, mean = -2, sd = 1)

Interactions of factors

Suppose fire alone reduces height by 2 centimeters \(\beta^{\text{fire}} = -2\) and drought alone reduces height by 4 centimeters \(\beta^{\text{drought}} = -4\). If a plant is exposed to both fire and drought, do you expect its height to be reduced by 6 centimeters ? Or more ? Or less ? Is both fire and drought different than the sum of the two effects ?

Maybe having both fire and drought is worse ? So the interaction between them is negative, let’s say \(\beta^{\text{fire},\text{drought}} = -2\).

Again, let’s write a small part our model for the tree heights in math:

\[\begin{align*} y_i &= \beta_0 \\ &\qquad + \beta^{\text{fire}} \text{Fire}_i + \beta^{\text{drought}} \text{Drought}_i + \beta^{\text{herbivory}} \text{Herbivory}_i \\ &\qquad + \beta^{\text{fire},\text{drought}} \text{Fire}_i \text{Drought}_i + \beta^{\text{fire},\text{herbivory}} \text{Fire}_i \text{Herbivory}_i + \beta^{\text{drought},\text{herbivory}} \text{Drought}_i \text{Herbivory}_i \\ &\qquad + \epsilon_i \end{align*}\]

Interactions can be visualized as lines that are not parallel. We plot only one possibile interaction from the model for simplicity:

Team activity (9 minutes): What does our model predict will be the height of a tree that gets fire and drought but not the herbivory treatment ?

An answer

This group gets \(\text{Fire}_i = 1\) and \(\text{Drought}_i = 1\), but \(\text{Herbivory}_i = 0\). So the model predicts their height will be:

\[\beta_0 + \beta^{\text{fire}}*1 + \beta^{\text{drought}}*1 + \beta^{\text{herbivory}}*0 + \beta^{\text{fire},\text{drought}} *1*1 + \beta^{\text{fire},\text{herbivory}} *1*0 + \beta^{\text{drought},\text{herbivory}} *1*0\] which is \(\beta_0 + \beta^{\text{fire}} + \beta^{\text{drought}} + \beta^{\text{fire},\text{drought}}\).

Team activity (11 minutes): Draw what you think is the interaction between fire and herbivory.

An answer

I was lazy and assume the other interactions were the same so I didn’t redo my drawing.

Our drawings are our best guesses for the interactions. But we don’t know them exactly, that is why we will collect data ! But maybe even before we collect data we know that the interactions are most likely negative. In other words, that having two treatments is even worse than the sum of both. We can express this prior information as:

\(\beta^{\text{fire},\text{drought}} \sim N(-2,1)\)

and same for the other interactions.

Let’s keep simulating to create fake data from our model’s story:

beta_fire_drought = rnorm(n = 1, mean = -2, sd = 1)
beta_fire_herbivory = rnorm(n = 1, mean = -2, sd = 1)
beta_drought_herbivory = rnorm(n = 1, mean = -2, sd = 1)

Species-specific effects

So far we’ve been talking in averages across 12 species of acacia trees:

Species_number Species
1 albhar
2 senbre
3 senmel
4 senpol
5 sensen
6 vacger
7 vacnil
8 vacrob
9 vacsey
10 vacsie
11 vactor
12 vacxan

Remember that we think the average height in the control group is most likely around 20 centimeters, and probably not less than 10 centimeters or more than 30 centimeters:

But our mathematical story might want to allow variation across species. Suppose average height is 20 centimeters, \(\beta_0 = 20\). Maybe senegalia brevispica is particularly short, so it is 5 centimeters below average, \(\beta^{\text{species}}_{\text{senegalia brevispica}} = -5\). Or maybe albizia harveyi is particularly tall, so it is 7 centimeters above average, \(\beta^{\text{species}}_{\text{albizia harveyi}} = 7\).

Maybe we think the shortest species is 8 centimeters below average and the tallest is 8 centimeters above average. Then we can express this mathematically with a Normal distribution with standard deviation 4 for the species-specific differences from the average:

\[\beta^{\text{species}}_s \sim N(0,4)\] In our mathematical story, each species’s height is the average height plus a choice from this Normal distribution:

Team activity (16 minutes): Use the R code below to generate some species heights. Do you get the same heights as your teammate ? Which species is the tallest ? Which is the shortest ?

beta_0 = 20
beta_species = rnorm(n = 12, mean = 0, sd = 4)
beta_0 + beta_species

An answer

set.seed(1)
beta_0 = 20
beta_species = rnorm(n = 12, mean = 0, sd = 4)
tibble(beta_0 + beta_species,
       Species = levels(as.factor(experimental_design$Species))) |>
  arrange(`beta_0 + beta_species`)
## # A tibble: 12 × 2
##    `beta_0 + beta_species` Species
##                      <dbl> <chr>  
##  1                    16.7 senmel 
##  2                    16.7 vacger 
##  3                    17.5 albhar 
##  4                    18.8 vacsie 
##  5                    20.7 senbre 
##  6                    21.3 sensen 
##  7                    21.6 vacxan 
##  8                    21.9 vacnil 
##  9                    22.3 vacsey 
## 10                    23.0 vacrob 
## 11                    26.0 vactor 
## 12                    26.4 senpol

It looks like in our simulated fake data world senpol (i.e. Senegalia polyacantha) is tallest and senmel (i.e. Senegalia melliferia) is shortest. This may not align with the real world because we have not incorporated species information into our model ! In the model, the species are exchangeable (p.5 of Bayesian Data Analysis).

Species-specific treatment effects

How much do the 12 species differ from each other in the effects of fire, drought, and herbivory ? Maybe about \(\pm 2\) centimeters ?

\[\beta^{\text{species,fire}}_s \sim N(0,1)\]

Suppose overall effect of fire is \(\beta^{\text{fire}} = -2\).

We next choose species-specific effects of fire \(\beta^{\text{fire}} + \beta^{\text{species,fire}}_s \sim N(\beta^{\text{fire}}, 1)\).


Multilevel modeling details and references

Maybe we don’t know how widely the species vary, so instead of choosing 4 and having \(\beta^{\text{species}}_s \sim N(0,4)\), we can add a level to our model to reflect our uncertainty:

\[\beta^{\text{species}}_s \sim N(0,\sigma_\text{species})\]

\[\sigma_\text{species} \sim \text{Exp}(\text{inverse-scale} = 1/4)\] We could have done this with the other parts of the model as well, but what makes this part special is that we have 12 species \(\beta^{\text{species}}_1,...,\beta^{\text{species}}_{12}\). So when we get data, we can learn about how these 12 species vary, informing \(\sigma_\text{species}\).

This technique is called multilevel modeling (also known as hierarchical modeling) and will be covered in a future course. If you want to read ahead, see:

sd_species = rexp(n = 1, rate = 1/4)
beta_species = rnorm(n = 12, mean = 0, sd = sd_species)
sd_species_fire = rexp(n = 1, rate = 1)
sd_species_drought = rexp(n = 1, rate = 1)
sd_species_herbivory = rexp(n = 1, rate = 1)
beta_species_fire = rnorm(n = 12, mean = 0, sd = sd_species_fire)
beta_species_drought = rnorm(n = 12, mean = 0, sd = sd_species_drought) 
beta_species_herbivory = rnorm(n = 12, mean = 0, sd = sd_species_herbivory)

Location effects

We do not directly aim to study the effect of plant location, but it might impact plant heights so we want to model it.

Team activity (8 minutes): How might a plant’s location in the experiment affect its growth ? Look back at the experimental design layout and think about each plant’s position.

An answer

If the soil or gradient changes across the experiment, we might expect the height in centimeters to change as the plant position moves south or east. Within a block, the position of the tree within the rainout shelter could affect its growth as well. Beyond the south/east directional effects, there could be block effects due to unintended differences in the rainout shelters.

Blocks 6-10 are located just downhill of Blocks 1-5. This might influence the soil in small but measurable ways.

Storms on the Nelson Mandela campus come from the east. This means that rain tends to be blown by the wind a bit. A small amount of rain is likely to blown into the east side of the experiment (Blocks 5 and 10). A tiny amount of rain might make it to the east sides of each Block between the roofs of the shelters.

The sun rises in the east, travels approximately north (the campus being south of the equator), and sets in the west. So the eastern side of the whole experiment (and potentially the individual blocks) might get more morning sun, and the west side might get more afternoon sun.

These small differences in sunlight and rain might influence our trees. We can attempt to measure and account for this.

We create variables \(\text{South}_i\), \(\text{East}_i\), \(\text{SouthInBlock}_i\), and \(\text{EastInBlock}_i\) for each tree \(i\) based on its plant position. Click the “Code” button below for the messy details.

SouthPlantsInSubplot = 3
EastPlantsInSubplot = 2
SouthSubplotsInBlock = 3
EastSubplotsInBlock = 3
SouthBlocks = 2
EastBlocks = 5

experimental_design = experimental_design %>% 
  mutate(
    # Calculate the position of the Plant in the Subplot:
    SouthPlantInSubplot_position = floor((PlantPosition-1)/EastPlantsInSubplot),
    EastPlantInSubplot_position = (PlantPosition-1) %% EastPlantsInSubplot,
    # Calculate the position of the Subplot in the Block:
    SouthSubplotInBlock_position = floor((Subplot-1)/EastSubplotsInBlock),
    EastSubplotInBlock_position = (Subplot-1) %% EastSubplotsInBlock,
    # Calculate the position of the Block in the whole experiment:
    SouthBlock_position = floor((Block-1)/EastBlocks),
    EastBlock_position = (Block-1) %% EastBlocks,
    
    # Combine the above to get the position of the Plant in Block:
    SouthPlantInBlock_position = SouthSubplotInBlock_position*SouthPlantsInSubplot + SouthPlantInSubplot_position,
    EastPlantInBlock_position = EastSubplotInBlock_position*EastPlantsInSubplot + EastPlantInSubplot_position,
    # Combine the above to get the position of the Plant in whole experiment:
    SouthPlant_position = SouthBlock_position*SouthSubplotsInBlock*SouthPlantsInSubplot + SouthPlantInBlock_position,
    EastPlant_position = EastBlock_position*EastSubplotsInBlock*EastPlantsInSubplot + EastPlantInBlock_position,
    
    # Center the position variables so they have mean zero:
    SouthInBlock = SouthPlantInBlock_position - mean(SouthPlantInBlock_position),
    EastInBlock = EastPlantInBlock_position - mean(EastPlantInBlock_position),
    South = SouthPlant_position - mean(SouthPlant_position),
    East = EastPlant_position - mean(EastPlant_position)
)

The model then allows these variables to affect height. For example, the coefficient of South position, \(\beta^{\text{South}}\), is the expected change in height in centimeters per change in South position. In general we expect these to be quite small, so our mathematical story gives them small standard deviations around zero:

beta_South = rnorm(n = 1, mean = 0, sd = 0.05)
beta_East = rnorm(n = 1, mean = 0, sd = 0.05)
beta_SouthInBlock = rnorm(n = 1, mean = 0, sd = 0.05)
beta_EastInBlock = rnorm(n = 1, mean = 0, sd = 0.05)

Beyond the South/East directional effects, there could be block effects due to unintended differences in the rainout shelters. Each block has a corn plant to capture environmental differences. Blocks with higher corn heights may have more hospitable conditions for acacia growth. For example, the coefficient of corn height, \(\beta^{\text{corn}}\), is the expected change in height in centimeters per change in corn height.

Team activity (9 minutes): For every centimeter increase in corn height, what do you expect to be the increase in acacia height for the acacias in that same block ? Express your guesses mathematically with a Normal distribution.

An answer

Corn grows more than acacia, so maybe 1 centimeter increase in corn corresponds to 0.5 centimeter advantage for acacia ? But I’m uncertain so I will put a standard deviation of 0.25 around it:

\(\beta^{\text{corn}} \sim N(0.5,0.25)\)

Again, we simulate to create fake data from our model’s story:

beta_corn = rnorm(n = 1, mean = 0.5, sd = 0.25)

Finally, there may still be a block effect we have not yet accounted for. We model this with a block effect, which we think is around \(\pm 4\) centimeters:

\[\beta^{\text{block}}_b \sim N(0,2)\]

sd_block = 2
beta_block = rnorm(n = 10, mean = 0, sd = sd_block)
Multilevel modeling details and references

As above for the species effects, we can add a level to our model to reflect our uncertainty:

\[\beta^{\text{block}}_b \sim N(0,\sigma_\text{block})\]

\[\sigma_\text{block} \sim \text{Exp}(\text{inverse-scale} = 1/2)\]

We have 10 blocks \(\beta^{\text{block}}_1,...,\beta^{\text{block}}_{10}\). So when we get data, we can learn about how these 10 blocks vary, informing \(\sigma_\text{block}\).

This technique is called multilevel modeling (also known as hierarchical modeling) and will be covered in a future course. If you want to read ahead, see:

Tree effects

Lastly, each tree may have its own specific effect, due to differences in individual seedlings. We model this with a tree effect \(\epsilon_i\), an “error” term that captures all the ways tree \(i\) is special that are not in our model.

Team activity (6 minutes): What might be special to a tree that isn’t in our model already ? In other words, what do you think is in the error term \(\epsilon_i\) ?

An answer

Each seedling will be genetically unique. Just like children in a family, they are all a little different ! These genetic differences might cause some trees to grow more quickly than others, resulting in differences in height.

Certain seedlings might be randomly eaten by insects or infected by various plant diseases. This might slow the growth of trees.

Some trees might get nibbled by rabbits or small antelopes (dikdik). This of course can change their height.

During the transplanting, some seedlings had their roots damaged. Other trees just were not happy to be transplanted, for unknown reasons. These transplant effects caused some individual seedlings to drop leaves or even cause branches to die. This can cause the height to decrease a bit.

Finally, measurement error ! Different scientists may get different results when measuring the same plant at the same time.

sd_tree = rhalft(n = 1, scale=2.5, nu=1) 
# brms default: 8. Parameters for specific families
# https://paul-buerkner.github.io/brms/reference/set_prior.html

Putting it together

We can now write our model for the tree heights in math (Greek letters !):

\[\begin{align*} y_i &= \beta_0 \\ &\qquad + \beta^{\text{fire}} \text{Fire}_i + \beta^{\text{drought}} \text{Drought}_i + \beta^{\text{herbivory}} \text{Herbivory}_i \\ &\qquad + \beta^{\text{fire},\text{drought}} \text{Fire}_i \text{Drought}_i + \beta^{\text{fire},\text{herbivory}} \text{Fire}_i \text{Herbivory}_i + \beta^{\text{drought},\text{herbivory}} \text{Drought}_i \text{Herbivory}_i \\ &\qquad + \beta^{\text{species}}_{s[i]} \\ &\qquad + \beta^{\text{species},\text{fire}}_{s[i]} \text{Fire}_i + \beta^{\text{species},\text{drought}}_{s[i]} \text{Drought}_i + \beta^{\text{species},\text{herbivory}}_{s[i]} \text{Herbivory}_i \\ &\qquad + \beta^{\text{SouthInBlock}} \text{SouthInBlock}_i + \beta^{\text{EastInBlock}} \text{EastInBlock}_i \\ &\qquad + \beta^{\text{South}} \text{South}_i + \beta^{\text{East}} \text{East}_i \\ &\qquad + \beta^{\text{corn}} \text{corn}_{b[i]} \\ &\qquad + \beta^{\text{block}}_{b[i]} \\ &\qquad + \epsilon_i \end{align*}\]

For now, let’s ignore the corn effect because we do not yet have corn height data.

Consider tree \(i = 1\) in the first row of data below:

tree Fire Drought Herbivory Species Species_number SouthInBlock EastInBlock South East Block is_corn
1 0 1 1 vacsey 9 -4.0102041 -2.5510204 -8.5102041 -14.5510204 1 0
2 0 1 1 senpol 4 -4.0102041 -1.5510204 -8.5102041 -13.5510204 1 0
3 0 1 1 vacger 6 -3.0102041 -2.5510204 -7.5102041 -14.5510204 1 0
4 0 1 1 albhar 1 -3.0102041 -1.5510204 -7.5102041 -13.5510204 1 0
5 0 1 1 senbre 2 -2.0102041 -2.5510204 -6.5102041 -14.5510204 1 0
6 0 1 1 sensen 5 -2.0102041 -1.5510204 -6.5102041 -13.5510204 1 0
7 0 0 1 albhar 1 -4.0102041 -0.5510204 -8.5102041 -12.5510204 1 0
8 0 0 1 vacrob 8 -4.0102041 0.4489796 -8.5102041 -11.5510204 1 0
9 0 0 1 senmel 3 -3.0102041 -0.5510204 -7.5102041 -12.5510204 1 0
10 0 0 1 vactor 11 -3.0102041 0.4489796 -7.5102041 -11.5510204 1 0
11 0 0 1 vacxan 12 -2.0102041 -0.5510204 -6.5102041 -12.5510204 1 0
12 0 0 1 vacnil 7 -2.0102041 0.4489796 -6.5102041 -11.5510204 1 0
13 0 0 0 albhar 1 -4.0102041 1.4489796 -8.5102041 -10.5510204 1 0
14 0 0 0 vacsie 10 -4.0102041 2.4489796 -8.5102041 -9.5510204 1 0
15 0 0 0 vacrob 8 -3.0102041 1.4489796 -7.5102041 -10.5510204 1 0
16 0 0 0 senmel 3 -3.0102041 2.4489796 -7.5102041 -9.5510204 1 0
17 0 0 0 senbre 2 -2.0102041 1.4489796 -6.5102041 -10.5510204 1 0
18 0 0 0 sensen 5 -2.0102041 2.4489796 -6.5102041 -9.5510204 1 0
19 0 1 0 vacxan 12 -1.0102041 -2.5510204 -5.5102041 -14.5510204 1 0
20 0 1 0 vactor 11 -1.0102041 -1.5510204 -5.5102041 -13.5510204 1 0
21 0 1 0 vacrob 8 -0.0102041 -2.5510204 -4.5102041 -14.5510204 1 0
22 0 1 0 senpol 4 -0.0102041 -1.5510204 -4.5102041 -13.5510204 1 0
23 0 1 0 senbre 2 0.9897959 -2.5510204 -3.5102041 -14.5510204 1 0
24 0 1 0 senmel 3 0.9897959 -1.5510204 -3.5102041 -13.5510204 1 0
25 0 0 1 sensen 5 -1.0102041 -0.5510204 -5.5102041 -12.5510204 1 0
26 0 0 1 senbre 2 -1.0102041 0.4489796 -5.5102041 -11.5510204 1 0
27 0 0 1 vacsie 10 -0.0102041 -0.5510204 -4.5102041 -12.5510204 1 0
28 0 0 1 vacsey 9 -0.0102041 0.4489796 -4.5102041 -11.5510204 1 0
29 0 0 1 senpol 4 0.9897959 -0.5510204 -3.5102041 -12.5510204 1 0
30 0 0 1 vacger 6 0.9897959 0.4489796 -3.5102041 -11.5510204 1 0
31 0 1 1 senmel 3 -1.0102041 1.4489796 -5.5102041 -10.5510204 1 0
32 0 1 1 vacrob 8 -1.0102041 2.4489796 -5.5102041 -9.5510204 1 0
33 0 1 1 vacnil 7 -0.0102041 1.4489796 -4.5102041 -10.5510204 1 0
34 0 1 1 vactor 11 -0.0102041 2.4489796 -4.5102041 -9.5510204 1 0
35 0 1 1 vacxan 12 0.9897959 1.4489796 -3.5102041 -10.5510204 1 0
36 0 1 1 vacsie 10 0.9897959 2.4489796 -3.5102041 -9.5510204 1 0
37 0 0 0 vacger 6 1.9897959 -2.5510204 -2.5102041 -14.5510204 1 1
38 0 0 0 vacnil 7 1.9897959 -0.5510204 -2.5102041 -12.5510204 1 0
39 0 0 0 senpol 4 1.9897959 0.4489796 -2.5102041 -11.5510204 1 0
40 0 0 0 vactor 11 2.9897959 -0.5510204 -1.5102041 -12.5510204 1 0
41 0 0 0 vacsey 9 2.9897959 0.4489796 -1.5102041 -11.5510204 1 0
42 0 0 0 vacger 6 3.9897959 -0.5510204 -0.5102041 -12.5510204 1 0
43 0 0 0 vacxan 12 3.9897959 0.4489796 -0.5102041 -11.5510204 1 0
44 0 1 0 sensen 5 1.9897959 1.4489796 -2.5102041 -10.5510204 1 0
45 0 1 0 vacsie 10 1.9897959 2.4489796 -2.5102041 -9.5510204 1 0
46 0 1 0 vacnil 7 2.9897959 1.4489796 -1.5102041 -10.5510204 1 0
47 0 1 0 vacger 6 2.9897959 2.4489796 -1.5102041 -9.5510204 1 0
48 0 1 0 albhar 1 3.9897959 1.4489796 -0.5102041 -10.5510204 1 0
49 0 1 0 vacsey 9 3.9897959 2.4489796 -0.5102041 -9.5510204 1 0
50 1 0 0 vacrob 8 -4.0102041 -2.5510204 -8.5102041 -8.5510204 2 1
51 1 0 0 vacxan 12 -4.0102041 -0.5510204 -8.5102041 -6.5510204 2 0
52 1 0 0 senbre 2 -4.0102041 0.4489796 -8.5102041 -5.5510204 2 0
53 1 0 0 sensen 5 -3.0102041 -0.5510204 -7.5102041 -6.5510204 2 0
54 1 0 0 vacrob 8 -3.0102041 0.4489796 -7.5102041 -5.5510204 2 0
55 1 0 0 senpol 4 -2.0102041 -0.5510204 -6.5102041 -6.5510204 2 0
56 1 0 0 vactor 11 -2.0102041 0.4489796 -6.5102041 -5.5510204 2 0
57 1 0 1 vacsey 9 -4.0102041 1.4489796 -8.5102041 -4.5510204 2 0
58 1 0 1 vacger 6 -4.0102041 2.4489796 -8.5102041 -3.5510204 2 0
59 1 0 1 senpol 4 -3.0102041 1.4489796 -7.5102041 -4.5510204 2 0
60 1 0 1 vactor 11 -3.0102041 2.4489796 -7.5102041 -3.5510204 2 0
61 1 0 1 sensen 5 -2.0102041 1.4489796 -6.5102041 -4.5510204 2 0
62 1 0 1 vacrob 8 -2.0102041 2.4489796 -6.5102041 -3.5510204 2 0
63 1 1 1 vacsie 10 -1.0102041 -2.5510204 -5.5102041 -8.5510204 2 0
64 1 1 1 vacger 6 -1.0102041 -1.5510204 -5.5102041 -7.5510204 2 0
65 1 1 1 senmel 3 -0.0102041 -2.5510204 -4.5102041 -8.5510204 2 0
66 1 1 1 vacnil 7 -0.0102041 -1.5510204 -4.5102041 -7.5510204 2 0
67 1 1 1 vacrob 8 0.9897959 -2.5510204 -3.5102041 -8.5510204 2 0
68 1 1 1 vacsey 9 0.9897959 -1.5510204 -3.5102041 -7.5510204 2 0
69 1 1 0 albhar 1 -1.0102041 -0.5510204 -5.5102041 -6.5510204 2 0
70 1 1 0 senmel 3 -1.0102041 0.4489796 -5.5102041 -5.5510204 2 0
71 1 1 0 vactor 11 -0.0102041 -0.5510204 -4.5102041 -6.5510204 2 0
72 1 1 0 vacrob 8 -0.0102041 0.4489796 -4.5102041 -5.5510204 2 0
73 1 1 0 vacnil 7 0.9897959 -0.5510204 -3.5102041 -6.5510204 2 0
74 1 1 0 vacger 6 0.9897959 0.4489796 -3.5102041 -5.5510204 2 0
75 1 1 0 senpol 4 -1.0102041 1.4489796 -5.5102041 -4.5510204 2 0
76 1 1 0 vacsie 10 -1.0102041 2.4489796 -5.5102041 -3.5510204 2 0
77 1 1 0 vacxan 12 -0.0102041 1.4489796 -4.5102041 -4.5510204 2 0
78 1 1 0 vacsey 9 -0.0102041 2.4489796 -4.5102041 -3.5510204 2 0
79 1 1 0 senbre 2 0.9897959 1.4489796 -3.5102041 -4.5510204 2 0
80 1 1 0 sensen 5 0.9897959 2.4489796 -3.5102041 -3.5510204 2 0
81 1 0 0 vacnil 7 1.9897959 -2.5510204 -2.5102041 -8.5510204 2 0
82 1 0 0 albhar 1 1.9897959 -1.5510204 -2.5102041 -7.5510204 2 0
83 1 0 0 vacsey 9 2.9897959 -2.5510204 -1.5102041 -8.5510204 2 0
84 1 0 0 vacger 6 2.9897959 -1.5510204 -1.5102041 -7.5510204 2 0
85 1 0 0 vacsie 10 3.9897959 -2.5510204 -0.5102041 -8.5510204 2 0
86 1 0 0 senmel 3 3.9897959 -1.5510204 -0.5102041 -7.5510204 2 0
87 1 0 1 vacsie 10 1.9897959 -0.5510204 -2.5102041 -6.5510204 2 0
88 1 0 1 vacxan 12 1.9897959 0.4489796 -2.5102041 -5.5510204 2 0
89 1 0 1 vacnil 7 2.9897959 -0.5510204 -1.5102041 -6.5510204 2 0
90 1 0 1 albhar 1 2.9897959 0.4489796 -1.5102041 -5.5510204 2 0
91 1 0 1 senmel 3 3.9897959 -0.5510204 -0.5102041 -6.5510204 2 0
92 1 0 1 senbre 2 3.9897959 0.4489796 -0.5102041 -5.5510204 2 0
93 1 1 1 senbre 2 1.9897959 1.4489796 -2.5102041 -4.5510204 2 0
94 1 1 1 sensen 5 1.9897959 2.4489796 -2.5102041 -3.5510204 2 0
95 1 1 1 vactor 11 2.9897959 1.4489796 -1.5102041 -4.5510204 2 0
96 1 1 1 senpol 4 2.9897959 2.4489796 -1.5102041 -3.5510204 2 0
97 1 1 1 vacxan 12 3.9897959 1.4489796 -0.5102041 -4.5510204 2 0
98 1 1 1 albhar 1 3.9897959 2.4489796 -0.5102041 -3.5510204 2 0
99 1 0 0 vacnil 7 -4.0102041 -2.5510204 -8.5102041 -2.5510204 3 0
100 1 0 0 vacsie 10 -4.0102041 -1.5510204 -8.5102041 -1.5510204 3 0
101 1 0 0 vacger 6 -3.0102041 -2.5510204 -7.5102041 -2.5510204 3 0
102 1 0 0 vactor 11 -3.0102041 -1.5510204 -7.5102041 -1.5510204 3 0
103 1 0 0 sensen 5 -2.0102041 -2.5510204 -6.5102041 -2.5510204 3 0
104 1 0 0 vacrob 8 -2.0102041 -1.5510204 -6.5102041 -1.5510204 3 0
105 1 0 1 albhar 1 -4.0102041 -0.5510204 -8.5102041 -0.5510204 3 0
106 1 0 1 vacnil 7 -4.0102041 0.4489796 -8.5102041 0.4489796 3 0
107 1 0 1 vactor 11 -3.0102041 -0.5510204 -7.5102041 -0.5510204 3 0
108 1 0 1 senpol 4 -3.0102041 0.4489796 -7.5102041 0.4489796 3 0
109 1 0 1 vacxan 12 -2.0102041 -0.5510204 -6.5102041 -0.5510204 3 0
110 1 0 1 sensen 5 -2.0102041 0.4489796 -6.5102041 0.4489796 3 0
111 1 1 1 vacsie 10 -4.0102041 1.4489796 -8.5102041 1.4489796 3 0
112 1 1 1 vacger 6 -4.0102041 2.4489796 -8.5102041 2.4489796 3 0
113 1 1 1 senmel 3 -3.0102041 1.4489796 -7.5102041 1.4489796 3 0
114 1 1 1 senbre 2 -3.0102041 2.4489796 -7.5102041 2.4489796 3 0
115 1 1 1 vacsey 9 -2.0102041 1.4489796 -6.5102041 1.4489796 3 0
116 1 1 1 albhar 1 -2.0102041 2.4489796 -6.5102041 2.4489796 3 0
117 1 0 0 senbre 2 -1.0102041 -2.5510204 -5.5102041 -2.5510204 3 1
118 1 1 0 senpol 4 -1.0102041 -0.5510204 -5.5102041 -0.5510204 3 0
119 1 1 0 albhar 1 -1.0102041 0.4489796 -5.5102041 0.4489796 3 0
120 1 1 0 senmel 3 -0.0102041 -0.5510204 -4.5102041 -0.5510204 3 0
121 1 1 0 senbre 2 -0.0102041 0.4489796 -4.5102041 0.4489796 3 0
122 1 1 0 vacsey 9 0.9897959 -0.5510204 -3.5102041 -0.5510204 3 0
123 1 1 0 vacnil 7 0.9897959 0.4489796 -3.5102041 0.4489796 3 0
124 1 1 1 vacxan 12 -1.0102041 1.4489796 -5.5102041 1.4489796 3 0
125 1 1 1 vacnil 7 -1.0102041 2.4489796 -5.5102041 2.4489796 3 0
126 1 1 1 vacrob 8 -0.0102041 1.4489796 -4.5102041 1.4489796 3 0
127 1 1 1 sensen 5 -0.0102041 2.4489796 -4.5102041 2.4489796 3 0
128 1 1 1 senpol 4 0.9897959 1.4489796 -3.5102041 1.4489796 3 0
129 1 1 1 vactor 11 0.9897959 2.4489796 -3.5102041 2.4489796 3 0
130 1 1 0 vacrob 8 1.9897959 -2.5510204 -2.5102041 -2.5510204 3 0
131 1 1 0 vacsie 10 1.9897959 -1.5510204 -2.5102041 -1.5510204 3 0
132 1 1 0 sensen 5 2.9897959 -2.5510204 -1.5102041 -2.5510204 3 0
133 1 1 0 vacger 6 2.9897959 -1.5510204 -1.5102041 -1.5510204 3 0
134 1 1 0 vacxan 12 3.9897959 -2.5510204 -0.5102041 -2.5510204 3 0
135 1 1 0 vactor 11 3.9897959 -1.5510204 -0.5102041 -1.5510204 3 0
136 1 0 1 vacger 6 1.9897959 -0.5510204 -2.5102041 -0.5510204 3 0
137 1 0 1 vacrob 8 1.9897959 0.4489796 -2.5102041 0.4489796 3 0
138 1 0 1 senmel 3 2.9897959 -0.5510204 -1.5102041 -0.5510204 3 0
139 1 0 1 senbre 2 2.9897959 0.4489796 -1.5102041 0.4489796 3 0
140 1 0 1 vacsie 10 3.9897959 -0.5510204 -0.5102041 -0.5510204 3 0
141 1 0 1 vacsey 9 3.9897959 0.4489796 -0.5102041 0.4489796 3 0
142 1 0 0 senmel 3 1.9897959 1.4489796 -2.5102041 1.4489796 3 0
143 1 0 0 vacxan 12 1.9897959 2.4489796 -2.5102041 2.4489796 3 0
144 1 0 0 albhar 1 2.9897959 1.4489796 -1.5102041 1.4489796 3 0
145 1 0 0 vacsey 9 2.9897959 2.4489796 -1.5102041 2.4489796 3 0
146 1 0 0 senpol 4 3.9897959 1.4489796 -0.5102041 1.4489796 3 0
147 1 0 0 senbre 2 3.9897959 2.4489796 -0.5102041 2.4489796 3 0
148 0 0 0 vacxan 12 -4.0102041 -2.5510204 -8.5102041 3.4489796 4 1
149 0 1 0 vacsey 9 -4.0102041 -0.5510204 -8.5102041 5.4489796 4 0
150 0 1 0 vacsie 10 -4.0102041 0.4489796 -8.5102041 6.4489796 4 0
151 0 1 0 vactor 11 -3.0102041 -0.5510204 -7.5102041 5.4489796 4 0
152 0 1 0 senbre 2 -3.0102041 0.4489796 -7.5102041 6.4489796 4 0
153 0 1 0 vacxan 12 -2.0102041 -0.5510204 -6.5102041 5.4489796 4 0
154 0 1 0 senmel 3 -2.0102041 0.4489796 -6.5102041 6.4489796 4 0
155 0 0 0 sensen 5 -4.0102041 1.4489796 -8.5102041 7.4489796 4 0
156 0 0 0 vacger 6 -4.0102041 2.4489796 -8.5102041 8.4489796 4 0
157 0 0 0 vacxan 12 -3.0102041 1.4489796 -7.5102041 7.4489796 4 0
158 0 0 0 vacnil 7 -3.0102041 2.4489796 -7.5102041 8.4489796 4 0
159 0 0 0 senpol 4 -2.0102041 1.4489796 -6.5102041 7.4489796 4 0
160 0 0 0 vactor 11 -2.0102041 2.4489796 -6.5102041 8.4489796 4 0
161 0 1 1 vacger 6 -1.0102041 -2.5510204 -5.5102041 3.4489796 4 0
162 0 1 1 senpol 4 -1.0102041 -1.5510204 -5.5102041 4.4489796 4 0
163 0 1 1 vacsey 9 -0.0102041 -2.5510204 -4.5102041 3.4489796 4 0
164 0 1 1 vacrob 8 -0.0102041 -1.5510204 -4.5102041 4.4489796 4 0
165 0 1 1 senmel 3 0.9897959 -2.5510204 -3.5102041 3.4489796 4 0
166 0 1 1 sensen 5 0.9897959 -1.5510204 -3.5102041 4.4489796 4 0
167 0 1 0 albhar 1 -1.0102041 -0.5510204 -5.5102041 5.4489796 4 0
168 0 1 0 vacrob 8 -1.0102041 0.4489796 -5.5102041 6.4489796 4 0
169 0 1 0 vacger 6 -0.0102041 -0.5510204 -4.5102041 5.4489796 4 0
170 0 1 0 senpol 4 -0.0102041 0.4489796 -4.5102041 6.4489796 4 0
171 0 1 0 vacnil 7 0.9897959 -0.5510204 -3.5102041 5.4489796 4 0
172 0 1 0 sensen 5 0.9897959 0.4489796 -3.5102041 6.4489796 4 0
173 0 0 1 senmel 3 -1.0102041 1.4489796 -5.5102041 7.4489796 4 0
174 0 0 1 vactor 11 -1.0102041 2.4489796 -5.5102041 8.4489796 4 0
175 0 0 1 vacrob 8 -0.0102041 1.4489796 -4.5102041 7.4489796 4 0
176 0 0 1 vacsey 9 -0.0102041 2.4489796 -4.5102041 8.4489796 4 0
177 0 0 1 senbre 2 0.9897959 1.4489796 -3.5102041 7.4489796 4 0
178 0 0 1 senpol 4 0.9897959 2.4489796 -3.5102041 8.4489796 4 0
179 0 0 1 vacsie 10 1.9897959 -2.5510204 -2.5102041 3.4489796 4 0
180 0 0 1 vacnil 7 1.9897959 -1.5510204 -2.5102041 4.4489796 4 0
181 0 0 1 vacxan 12 2.9897959 -2.5510204 -1.5102041 3.4489796 4 0
182 0 0 1 albhar 1 2.9897959 -1.5510204 -1.5102041 4.4489796 4 0
183 0 0 1 vacger 6 3.9897959 -2.5510204 -0.5102041 3.4489796 4 0
184 0 0 1 sensen 5 3.9897959 -1.5510204 -0.5102041 4.4489796 4 0
185 0 0 0 senbre 2 1.9897959 -0.5510204 -2.5102041 5.4489796 4 0
186 0 0 0 albhar 1 1.9897959 0.4489796 -2.5102041 6.4489796 4 0
187 0 0 0 vacrob 8 2.9897959 -0.5510204 -1.5102041 5.4489796 4 0
188 0 0 0 senmel 3 2.9897959 0.4489796 -1.5102041 6.4489796 4 0
189 0 0 0 vacsey 9 3.9897959 -0.5510204 -0.5102041 5.4489796 4 0
190 0 0 0 vacsie 10 3.9897959 0.4489796 -0.5102041 6.4489796 4 0
191 0 1 1 vacsie 10 1.9897959 1.4489796 -2.5102041 7.4489796 4 0
192 0 1 1 vactor 11 1.9897959 2.4489796 -2.5102041 8.4489796 4 0
193 0 1 1 senbre 2 2.9897959 1.4489796 -1.5102041 7.4489796 4 0
194 0 1 1 albhar 1 2.9897959 2.4489796 -1.5102041 8.4489796 4 0
195 0 1 1 vacxan 12 3.9897959 1.4489796 -0.5102041 7.4489796 4 0
196 0 1 1 vacnil 7 3.9897959 2.4489796 -0.5102041 8.4489796 4 0
197 0 0 0 vacsey 9 -4.0102041 -2.5510204 -8.5102041 9.4489796 5 0
198 0 0 0 vacxan 12 -4.0102041 -1.5510204 -8.5102041 10.4489796 5 0
199 0 0 0 sensen 5 -3.0102041 -2.5510204 -7.5102041 9.4489796 5 0
200 0 0 0 senmel 3 -3.0102041 -1.5510204 -7.5102041 10.4489796 5 0
201 0 0 0 vacnil 7 -2.0102041 -2.5510204 -6.5102041 9.4489796 5 0
202 0 0 0 vactor 11 -2.0102041 -1.5510204 -6.5102041 10.4489796 5 0
203 0 0 1 vacnil 7 -4.0102041 -0.5510204 -8.5102041 11.4489796 5 0
204 0 0 1 senbre 2 -4.0102041 0.4489796 -8.5102041 12.4489796 5 0
205 0 0 1 vactor 11 -3.0102041 -0.5510204 -7.5102041 11.4489796 5 0
206 0 0 1 vacrob 8 -3.0102041 0.4489796 -7.5102041 12.4489796 5 0
207 0 0 1 vacxan 12 -2.0102041 -0.5510204 -6.5102041 11.4489796 5 0
208 0 0 1 vacsey 9 -2.0102041 0.4489796 -6.5102041 12.4489796 5 0
209 0 1 0 vacrob 8 -4.0102041 1.4489796 -8.5102041 13.4489796 5 0
210 0 1 0 vactor 11 -4.0102041 2.4489796 -8.5102041 14.4489796 5 0
211 0 1 0 vacsey 9 -3.0102041 1.4489796 -7.5102041 13.4489796 5 0
212 0 1 0 vacsie 10 -3.0102041 2.4489796 -7.5102041 14.4489796 5 0
213 0 1 0 senmel 3 -2.0102041 1.4489796 -6.5102041 13.4489796 5 0
214 0 1 0 vacger 6 -2.0102041 2.4489796 -6.5102041 14.4489796 5 0
215 0 0 0 vacrob 8 -1.0102041 -2.5510204 -5.5102041 9.4489796 5 0
216 0 0 0 vacger 6 -1.0102041 -1.5510204 -5.5102041 10.4489796 5 0
217 0 0 0 senpol 4 -0.0102041 -2.5510204 -4.5102041 9.4489796 5 0
218 0 0 0 vacsie 10 -0.0102041 -1.5510204 -4.5102041 10.4489796 5 0
219 0 0 0 albhar 1 0.9897959 -2.5510204 -3.5102041 9.4489796 5 0
220 0 0 0 senbre 2 0.9897959 -1.5510204 -3.5102041 10.4489796 5 0
221 0 1 1 vacnil 7 -1.0102041 -0.5510204 -5.5102041 11.4489796 5 0
222 0 1 1 albhar 1 -1.0102041 0.4489796 -5.5102041 12.4489796 5 0
223 0 1 1 senpol 4 -0.0102041 -0.5510204 -4.5102041 11.4489796 5 0
224 0 1 1 sensen 5 -0.0102041 0.4489796 -4.5102041 12.4489796 5 0
225 0 1 1 vactor 11 0.9897959 -0.5510204 -3.5102041 11.4489796 5 0
226 0 1 1 senmel 3 0.9897959 0.4489796 -3.5102041 12.4489796 5 0
227 0 0 1 senpol 4 -1.0102041 1.4489796 -5.5102041 13.4489796 5 0
228 0 0 1 vacsie 10 -1.0102041 2.4489796 -5.5102041 14.4489796 5 0
229 0 0 1 vacger 6 -0.0102041 1.4489796 -4.5102041 13.4489796 5 0
230 0 0 1 albhar 1 -0.0102041 2.4489796 -4.5102041 14.4489796 5 0
231 0 0 1 senmel 3 0.9897959 1.4489796 -3.5102041 13.4489796 5 0
232 0 0 1 sensen 5 0.9897959 2.4489796 -3.5102041 14.4489796 5 0
233 0 1 1 vacsie 10 1.9897959 -2.5510204 -2.5102041 9.4489796 5 0
234 0 1 1 vacsey 9 1.9897959 -1.5510204 -2.5102041 10.4489796 5 0
235 0 1 1 vacxan 12 2.9897959 -2.5510204 -1.5102041 9.4489796 5 0
236 0 1 1 vacrob 8 2.9897959 -1.5510204 -1.5102041 10.4489796 5 0
237 0 1 1 vacger 6 3.9897959 -2.5510204 -0.5102041 9.4489796 5 0
238 0 1 1 senbre 2 3.9897959 -1.5510204 -0.5102041 10.4489796 5 0
239 0 0 0 vactor 11 1.9897959 -0.5510204 -2.5102041 11.4489796 5 1
240 0 1 0 sensen 5 1.9897959 1.4489796 -2.5102041 13.4489796 5 0
241 0 1 0 albhar 1 1.9897959 2.4489796 -2.5102041 14.4489796 5 0
242 0 1 0 vacxan 12 2.9897959 1.4489796 -1.5102041 13.4489796 5 0
243 0 1 0 senpol 4 2.9897959 2.4489796 -1.5102041 14.4489796 5 0
244 0 1 0 senbre 2 3.9897959 1.4489796 -0.5102041 13.4489796 5 0
245 0 1 0 vacnil 7 3.9897959 2.4489796 -0.5102041 14.4489796 5 0
246 0 0 0 sensen 5 -4.0102041 -2.5510204 0.4897959 -14.5510204 6 0
247 0 0 0 albhar 1 -4.0102041 -1.5510204 0.4897959 -13.5510204 6 0
248 0 0 0 senbre 2 -3.0102041 -2.5510204 1.4897959 -14.5510204 6 0
249 0 0 0 vactor 11 -3.0102041 -1.5510204 1.4897959 -13.5510204 6 0
250 0 0 0 vacrob 8 -2.0102041 -2.5510204 2.4897959 -14.5510204 6 0
251 0 0 0 vacxan 12 -2.0102041 -1.5510204 2.4897959 -13.5510204 6 0
252 0 0 1 vacxan 12 -4.0102041 -0.5510204 0.4897959 -12.5510204 6 0
253 0 0 1 vacsie 10 -4.0102041 0.4489796 0.4897959 -11.5510204 6 0
254 0 0 1 senbre 2 -3.0102041 -0.5510204 1.4897959 -12.5510204 6 0
255 0 0 1 albhar 1 -3.0102041 0.4489796 1.4897959 -11.5510204 6 0
256 0 0 1 vacsey 9 -2.0102041 -0.5510204 2.4897959 -12.5510204 6 0
257 0 0 1 senpol 4 -2.0102041 0.4489796 2.4897959 -11.5510204 6 0
258 0 1 1 vacsey 9 -4.0102041 1.4489796 0.4897959 -10.5510204 6 0
259 0 1 1 vacnil 7 -4.0102041 2.4489796 0.4897959 -9.5510204 6 0
260 0 1 1 vactor 11 -3.0102041 1.4489796 1.4897959 -10.5510204 6 0
261 0 1 1 senpol 4 -3.0102041 2.4489796 1.4897959 -9.5510204 6 0
262 0 1 1 senmel 3 -2.0102041 1.4489796 2.4897959 -10.5510204 6 0
263 0 1 1 sensen 5 -2.0102041 2.4489796 2.4897959 -9.5510204 6 0
264 0 1 0 sensen 5 -1.0102041 -2.5510204 3.4897959 -14.5510204 6 0
265 0 1 0 vacger 6 -1.0102041 -1.5510204 3.4897959 -13.5510204 6 0
266 0 1 0 senpol 4 -0.0102041 -2.5510204 4.4897959 -14.5510204 6 0
267 0 1 0 vacsie 10 -0.0102041 -1.5510204 4.4897959 -13.5510204 6 0
268 0 1 0 vacxan 12 0.9897959 -2.5510204 5.4897959 -14.5510204 6 0
269 0 1 0 vactor 11 0.9897959 -1.5510204 5.4897959 -13.5510204 6 0
270 0 0 0 senmel 3 -1.0102041 -0.5510204 3.4897959 -12.5510204 6 1
271 0 0 1 vacger 6 -1.0102041 1.4489796 3.4897959 -10.5510204 6 0
272 0 0 1 sensen 5 -1.0102041 2.4489796 3.4897959 -9.5510204 6 0
273 0 0 1 vacrob 8 -0.0102041 1.4489796 4.4897959 -10.5510204 6 0
274 0 0 1 vactor 11 -0.0102041 2.4489796 4.4897959 -9.5510204 6 0
275 0 0 1 vacnil 7 0.9897959 1.4489796 5.4897959 -10.5510204 6 0
276 0 0 1 senmel 3 0.9897959 2.4489796 5.4897959 -9.5510204 6 0
277 0 1 0 senbre 2 1.9897959 -2.5510204 6.4897959 -14.5510204 6 0
278 0 1 0 vacsey 9 1.9897959 -1.5510204 6.4897959 -13.5510204 6 0
279 0 1 0 vacrob 8 2.9897959 -2.5510204 7.4897959 -14.5510204 6 0
280 0 1 0 albhar 1 2.9897959 -1.5510204 7.4897959 -13.5510204 6 0
281 0 1 0 vacnil 7 3.9897959 -2.5510204 8.4897959 -14.5510204 6 0
282 0 1 0 senmel 3 3.9897959 -1.5510204 8.4897959 -13.5510204 6 0
283 0 1 1 vacrob 8 1.9897959 -0.5510204 6.4897959 -12.5510204 6 0
284 0 1 1 vacger 6 1.9897959 0.4489796 6.4897959 -11.5510204 6 0
285 0 1 1 albhar 1 2.9897959 -0.5510204 7.4897959 -12.5510204 6 0
286 0 1 1 vacxan 12 2.9897959 0.4489796 7.4897959 -11.5510204 6 0
287 0 1 1 vacsie 10 3.9897959 -0.5510204 8.4897959 -12.5510204 6 0
288 0 1 1 senbre 2 3.9897959 0.4489796 8.4897959 -11.5510204 6 0
289 0 0 0 vacsie 10 1.9897959 1.4489796 6.4897959 -10.5510204 6 0
290 0 0 0 senmel 3 1.9897959 2.4489796 6.4897959 -9.5510204 6 0
291 0 0 0 vacger 6 2.9897959 1.4489796 7.4897959 -10.5510204 6 0
292 0 0 0 vacnil 7 2.9897959 2.4489796 7.4897959 -9.5510204 6 0
293 0 0 0 senpol 4 3.9897959 1.4489796 8.4897959 -10.5510204 6 0
294 0 0 0 vacsey 9 3.9897959 2.4489796 8.4897959 -9.5510204 6 0
295 1 1 0 vacsey 9 -4.0102041 -2.5510204 0.4897959 -8.5510204 7 0
296 1 1 0 senmel 3 -4.0102041 -1.5510204 0.4897959 -7.5510204 7 0
297 1 1 0 vacrob 8 -3.0102041 -2.5510204 1.4897959 -8.5510204 7 0
298 1 1 0 vacger 6 -3.0102041 -1.5510204 1.4897959 -7.5510204 7 0
299 1 1 0 senpol 4 -2.0102041 -2.5510204 2.4897959 -8.5510204 7 0
300 1 1 0 vacnil 7 -2.0102041 -1.5510204 2.4897959 -7.5510204 7 0
301 1 0 0 senmel 3 -4.0102041 -0.5510204 0.4897959 -6.5510204 7 0
302 1 0 0 vacsie 10 -4.0102041 0.4489796 0.4897959 -5.5510204 7 0
303 1 0 0 vacnil 7 -3.0102041 -0.5510204 1.4897959 -6.5510204 7 0
304 1 0 0 vacger 6 -3.0102041 0.4489796 1.4897959 -5.5510204 7 0
305 1 0 0 albhar 1 -2.0102041 -0.5510204 2.4897959 -6.5510204 7 0
306 1 0 0 vacrob 8 -2.0102041 0.4489796 2.4897959 -5.5510204 7 0
307 1 0 1 albhar 1 -4.0102041 1.4489796 0.4897959 -4.5510204 7 0
308 1 0 1 vactor 11 -4.0102041 2.4489796 0.4897959 -3.5510204 7 0
309 1 0 1 vacsey 9 -3.0102041 1.4489796 1.4897959 -4.5510204 7 0
310 1 0 1 vacger 6 -3.0102041 2.4489796 1.4897959 -3.5510204 7 0
311 1 0 1 vacsie 10 -2.0102041 1.4489796 2.4897959 -4.5510204 7 0
312 1 0 1 vacnil 7 -2.0102041 2.4489796 2.4897959 -3.5510204 7 0
313 1 0 0 vactor 11 -1.0102041 -2.5510204 3.4897959 -8.5510204 7 0
314 1 0 0 vacsey 9 -1.0102041 -1.5510204 3.4897959 -7.5510204 7 0
315 1 0 0 vacxan 12 -0.0102041 -2.5510204 4.4897959 -8.5510204 7 0
316 1 0 0 senpol 4 -0.0102041 -1.5510204 4.4897959 -7.5510204 7 0
317 1 0 0 sensen 5 0.9897959 -2.5510204 5.4897959 -8.5510204 7 0
318 1 0 0 senbre 2 0.9897959 -1.5510204 5.4897959 -7.5510204 7 0
319 1 1 1 senpol 4 -1.0102041 -0.5510204 3.4897959 -6.5510204 7 0
320 1 1 1 senbre 2 -1.0102041 0.4489796 3.4897959 -5.5510204 7 0
321 1 1 1 sensen 5 -0.0102041 -0.5510204 4.4897959 -6.5510204 7 0
322 1 1 1 vactor 11 -0.0102041 0.4489796 4.4897959 -5.5510204 7 0
323 1 1 1 vacxan 12 0.9897959 -0.5510204 5.4897959 -6.5510204 7 0
324 1 1 1 senmel 3 0.9897959 0.4489796 5.4897959 -5.5510204 7 0
325 1 1 1 vacger 6 -1.0102041 1.4489796 3.4897959 -4.5510204 7 0
326 1 1 1 albhar 1 -1.0102041 2.4489796 3.4897959 -3.5510204 7 0
327 1 1 1 vacnil 7 -0.0102041 1.4489796 4.4897959 -4.5510204 7 0
328 1 1 1 vacrob 8 -0.0102041 2.4489796 4.4897959 -3.5510204 7 0
329 1 1 1 vacsie 10 0.9897959 1.4489796 5.4897959 -4.5510204 7 0
330 1 1 1 vacsey 9 0.9897959 2.4489796 5.4897959 -3.5510204 7 0
331 1 0 1 sensen 5 1.9897959 -2.5510204 6.4897959 -8.5510204 7 0
332 1 0 1 senpol 4 1.9897959 -1.5510204 6.4897959 -7.5510204 7 0
333 1 0 1 vacrob 8 2.9897959 -2.5510204 7.4897959 -8.5510204 7 0
334 1 0 1 senbre 2 2.9897959 -1.5510204 7.4897959 -7.5510204 7 0
335 1 0 1 vacxan 12 3.9897959 -2.5510204 8.4897959 -8.5510204 7 0
336 1 0 1 senmel 3 3.9897959 -1.5510204 8.4897959 -7.5510204 7 0
337 1 0 0 albhar 1 1.9897959 -0.5510204 6.4897959 -6.5510204 7 1
338 1 1 0 sensen 5 1.9897959 1.4489796 6.4897959 -4.5510204 7 0
339 1 1 0 vactor 11 1.9897959 2.4489796 6.4897959 -3.5510204 7 0
340 1 1 0 albhar 1 2.9897959 1.4489796 7.4897959 -4.5510204 7 0
341 1 1 0 vacsie 10 2.9897959 2.4489796 7.4897959 -3.5510204 7 0
342 1 1 0 senbre 2 3.9897959 1.4489796 8.4897959 -4.5510204 7 0
343 1 1 0 vacxan 12 3.9897959 2.4489796 8.4897959 -3.5510204 7 0
344 1 1 0 vacger 6 -4.0102041 -2.5510204 0.4897959 -2.5510204 8 0
345 1 1 0 albhar 1 -4.0102041 -1.5510204 0.4897959 -1.5510204 8 0
346 1 1 0 senpol 4 -3.0102041 -2.5510204 1.4897959 -2.5510204 8 0
347 1 1 0 vacsie 10 -3.0102041 -1.5510204 1.4897959 -1.5510204 8 0
348 1 1 0 senbre 2 -2.0102041 -2.5510204 2.4897959 -2.5510204 8 0
349 1 1 0 sensen 5 -2.0102041 -1.5510204 2.4897959 -1.5510204 8 0
350 1 0 0 sensen 5 -4.0102041 -0.5510204 0.4897959 -0.5510204 8 0
351 1 0 0 senmel 3 -4.0102041 0.4489796 0.4897959 0.4489796 8 0
352 1 0 0 senpol 4 -3.0102041 -0.5510204 1.4897959 -0.5510204 8 0
353 1 0 0 vacsey 9 -3.0102041 0.4489796 1.4897959 0.4489796 8 0
354 1 0 0 vacger 6 -2.0102041 -0.5510204 2.4897959 -0.5510204 8 0
355 1 0 0 vactor 11 -2.0102041 0.4489796 2.4897959 0.4489796 8 0
356 1 0 0 vacsey 9 -4.0102041 1.4489796 0.4897959 1.4489796 8 1
357 1 1 0 vacsey 9 -1.0102041 -2.5510204 3.4897959 -2.5510204 8 0
358 1 1 0 vactor 11 -1.0102041 -1.5510204 3.4897959 -1.5510204 8 0
359 1 1 0 vacnil 7 -0.0102041 -2.5510204 4.4897959 -2.5510204 8 0
360 1 1 0 vacrob 8 -0.0102041 -1.5510204 4.4897959 -1.5510204 8 0
361 1 1 0 senmel 3 0.9897959 -2.5510204 5.4897959 -2.5510204 8 0
362 1 1 0 vacxan 12 0.9897959 -1.5510204 5.4897959 -1.5510204 8 0
363 1 0 1 senbre 2 -1.0102041 -0.5510204 3.4897959 -0.5510204 8 0
364 1 0 1 senpol 4 -1.0102041 0.4489796 3.4897959 0.4489796 8 0
365 1 0 1 senmel 3 -0.0102041 -0.5510204 4.4897959 -0.5510204 8 0
366 1 0 1 vactor 11 -0.0102041 0.4489796 4.4897959 0.4489796 8 0
367 1 0 1 vacger 6 0.9897959 -0.5510204 5.4897959 -0.5510204 8 0
368 1 0 1 vacsie 10 0.9897959 0.4489796 5.4897959 0.4489796 8 0
369 1 1 1 senbre 2 -1.0102041 1.4489796 3.4897959 1.4489796 8 0
370 1 1 1 sensen 5 -1.0102041 2.4489796 3.4897959 2.4489796 8 0
371 1 1 1 vacsie 10 -0.0102041 1.4489796 4.4897959 1.4489796 8 0
372 1 1 1 vacrob 8 -0.0102041 2.4489796 4.4897959 2.4489796 8 0
373 1 1 1 vacxan 12 0.9897959 1.4489796 5.4897959 1.4489796 8 0
374 1 1 1 senpol 4 0.9897959 2.4489796 5.4897959 2.4489796 8 0
375 1 0 1 vacxan 12 1.9897959 -2.5510204 6.4897959 -2.5510204 8 0
376 1 0 1 albhar 1 1.9897959 -1.5510204 6.4897959 -1.5510204 8 0
377 1 0 1 vacrob 8 2.9897959 -2.5510204 7.4897959 -2.5510204 8 0
378 1 0 1 sensen 5 2.9897959 -1.5510204 7.4897959 -1.5510204 8 0
379 1 0 1 vacnil 7 3.9897959 -2.5510204 8.4897959 -2.5510204 8 0
380 1 0 1 vacsey 9 3.9897959 -1.5510204 8.4897959 -1.5510204 8 0
381 1 1 1 albhar 1 1.9897959 -0.5510204 6.4897959 -0.5510204 8 0
382 1 1 1 vacnil 7 1.9897959 0.4489796 6.4897959 0.4489796 8 0
383 1 1 1 vacger 6 2.9897959 -0.5510204 7.4897959 -0.5510204 8 0
384 1 1 1 vactor 11 2.9897959 0.4489796 7.4897959 0.4489796 8 0
385 1 1 1 senmel 3 3.9897959 -0.5510204 8.4897959 -0.5510204 8 0
386 1 1 1 vacsey 9 3.9897959 0.4489796 8.4897959 0.4489796 8 0
387 1 0 0 vacrob 8 1.9897959 1.4489796 6.4897959 1.4489796 8 0
388 1 0 0 vacsie 10 1.9897959 2.4489796 6.4897959 2.4489796 8 0
389 1 0 0 senbre 2 2.9897959 1.4489796 7.4897959 1.4489796 8 0
390 1 0 0 albhar 1 2.9897959 2.4489796 7.4897959 2.4489796 8 0
391 1 0 0 vacnil 7 3.9897959 1.4489796 8.4897959 1.4489796 8 0
392 1 0 0 vacxan 12 3.9897959 2.4489796 8.4897959 2.4489796 8 0
393 0 1 0 senmel 3 -4.0102041 -2.5510204 0.4897959 3.4489796 9 0
394 0 1 0 vacxan 12 -4.0102041 -1.5510204 0.4897959 4.4489796 9 0
395 0 1 0 albhar 1 -3.0102041 -2.5510204 1.4897959 3.4489796 9 0
396 0 1 0 senbre 2 -3.0102041 -1.5510204 1.4897959 4.4489796 9 0
397 0 1 0 vacrob 8 -2.0102041 -2.5510204 2.4897959 3.4489796 9 0
398 0 1 0 vacnil 7 -2.0102041 -1.5510204 2.4897959 4.4489796 9 0
399 0 0 1 vacsey 9 -4.0102041 -0.5510204 0.4897959 5.4489796 9 0
400 0 0 1 sensen 5 -4.0102041 0.4489796 0.4897959 6.4489796 9 0
401 0 0 1 albhar 1 -3.0102041 -0.5510204 1.4897959 5.4489796 9 0
402 0 0 1 senpol 4 -3.0102041 0.4489796 1.4897959 6.4489796 9 0
403 0 0 1 senmel 3 -2.0102041 -0.5510204 2.4897959 5.4489796 9 0
404 0 0 1 vacger 6 -2.0102041 0.4489796 2.4897959 6.4489796 9 0
405 0 0 0 sensen 5 -4.0102041 1.4489796 0.4897959 7.4489796 9 1
406 0 0 0 sensen 5 -1.0102041 -2.5510204 3.4897959 3.4489796 9 0
407 0 0 0 vacsie 10 -1.0102041 -1.5510204 3.4897959 4.4489796 9 0
408 0 0 0 vacrob 8 -0.0102041 -2.5510204 4.4897959 3.4489796 9 0
409 0 0 0 senbre 2 -0.0102041 -1.5510204 4.4897959 4.4489796 9 0
410 0 0 0 albhar 1 0.9897959 -2.5510204 5.4897959 3.4489796 9 0
411 0 0 0 vactor 11 0.9897959 -1.5510204 5.4897959 4.4489796 9 0
412 0 1 1 vacxan 12 -1.0102041 -0.5510204 3.4897959 5.4489796 9 0
413 0 1 1 albhar 1 -1.0102041 0.4489796 3.4897959 6.4489796 9 0
414 0 1 1 sensen 5 -0.0102041 -0.5510204 4.4897959 5.4489796 9 0
415 0 1 1 vacnil 7 -0.0102041 0.4489796 4.4897959 6.4489796 9 0
416 0 1 1 vacsey 9 0.9897959 -0.5510204 5.4897959 5.4489796 9 0
417 0 1 1 vactor 11 0.9897959 0.4489796 5.4897959 6.4489796 9 0
418 0 1 1 senmel 3 -1.0102041 1.4489796 3.4897959 7.4489796 9 0
419 0 1 1 senbre 2 -1.0102041 2.4489796 3.4897959 8.4489796 9 0
420 0 1 1 vacger 6 -0.0102041 1.4489796 4.4897959 7.4489796 9 0
421 0 1 1 senpol 4 -0.0102041 2.4489796 4.4897959 8.4489796 9 0
422 0 1 1 vacrob 8 0.9897959 1.4489796 5.4897959 7.4489796 9 0
423 0 1 1 vacsie 10 0.9897959 2.4489796 5.4897959 8.4489796 9 0
424 0 0 0 vacger 6 1.9897959 -2.5510204 6.4897959 3.4489796 9 0
425 0 0 0 senpol 4 1.9897959 -1.5510204 6.4897959 4.4489796 9 0
426 0 0 0 vacsey 9 2.9897959 -2.5510204 7.4897959 3.4489796 9 0
427 0 0 0 vacxan 12 2.9897959 -1.5510204 7.4897959 4.4489796 9 0
428 0 0 0 senmel 3 3.9897959 -2.5510204 8.4897959 3.4489796 9 0
429 0 0 0 vacnil 7 3.9897959 -1.5510204 8.4897959 4.4489796 9 0
430 0 1 0 senpol 4 1.9897959 -0.5510204 6.4897959 5.4489796 9 0
431 0 1 0 vacsey 9 1.9897959 0.4489796 6.4897959 6.4489796 9 0
432 0 1 0 vactor 11 2.9897959 -0.5510204 7.4897959 5.4489796 9 0
433 0 1 0 sensen 5 2.9897959 0.4489796 7.4897959 6.4489796 9 0
434 0 1 0 vacsie 10 3.9897959 -0.5510204 8.4897959 5.4489796 9 0
435 0 1 0 vacger 6 3.9897959 0.4489796 8.4897959 6.4489796 9 0
436 0 0 1 vactor 11 1.9897959 1.4489796 6.4897959 7.4489796 9 0
437 0 0 1 vacxan 12 1.9897959 2.4489796 6.4897959 8.4489796 9 0
438 0 0 1 vacsie 10 2.9897959 1.4489796 7.4897959 7.4489796 9 0
439 0 0 1 senbre 2 2.9897959 2.4489796 7.4897959 8.4489796 9 0
440 0 0 1 vacrob 8 3.9897959 1.4489796 8.4897959 7.4489796 9 0
441 0 0 1 vacnil 7 3.9897959 2.4489796 8.4897959 8.4489796 9 0
442 1 1 0 senmel 3 -4.0102041 -2.5510204 0.4897959 9.4489796 10 0
443 1 1 0 vactor 11 -4.0102041 -1.5510204 0.4897959 10.4489796 10 0
444 1 1 0 albhar 1 -3.0102041 -2.5510204 1.4897959 9.4489796 10 0
445 1 1 0 vacsie 10 -3.0102041 -1.5510204 1.4897959 10.4489796 10 0
446 1 1 0 vacxan 12 -2.0102041 -2.5510204 2.4897959 9.4489796 10 0
447 1 1 0 vacger 6 -2.0102041 -1.5510204 2.4897959 10.4489796 10 0
448 1 0 0 senpol 4 -4.0102041 -0.5510204 0.4897959 11.4489796 10 0
449 1 0 0 senmel 3 -4.0102041 0.4489796 0.4897959 12.4489796 10 0
450 1 0 0 sensen 5 -3.0102041 -0.5510204 1.4897959 11.4489796 10 0
451 1 0 0 senbre 2 -3.0102041 0.4489796 1.4897959 12.4489796 10 0
452 1 0 0 vacsie 10 -2.0102041 -0.5510204 2.4897959 11.4489796 10 0
453 1 0 0 vacrob 8 -2.0102041 0.4489796 2.4897959 12.4489796 10 0
454 1 0 1 vacrob 8 -4.0102041 1.4489796 0.4897959 13.4489796 10 0
455 1 0 1 albhar 1 -4.0102041 2.4489796 0.4897959 14.4489796 10 0
456 1 0 1 senmel 3 -3.0102041 1.4489796 1.4897959 13.4489796 10 0
457 1 0 1 vactor 11 -3.0102041 2.4489796 1.4897959 14.4489796 10 0
458 1 0 1 vacsie 10 -2.0102041 1.4489796 2.4897959 13.4489796 10 0
459 1 0 1 vacnil 7 -2.0102041 2.4489796 2.4897959 14.4489796 10 0
460 1 0 0 vacnil 7 -1.0102041 -2.5510204 3.4897959 9.4489796 10 1
461 1 0 1 vacger 6 -1.0102041 -0.5510204 3.4897959 11.4489796 10 0
462 1 0 1 vacxan 12 -1.0102041 0.4489796 3.4897959 12.4489796 10 0
463 1 0 1 senpol 4 -0.0102041 -0.5510204 4.4897959 11.4489796 10 0
464 1 0 1 sensen 5 -0.0102041 0.4489796 4.4897959 12.4489796 10 0
465 1 0 1 vacsey 9 0.9897959 -0.5510204 5.4897959 11.4489796 10 0
466 1 0 1 senbre 2 0.9897959 0.4489796 5.4897959 12.4489796 10 0
467 1 1 0 sensen 5 -1.0102041 1.4489796 3.4897959 13.4489796 10 0
468 1 1 0 vacrob 8 -1.0102041 2.4489796 3.4897959 14.4489796 10 0
469 1 1 0 vacsey 9 -0.0102041 1.4489796 4.4897959 13.4489796 10 0
470 1 1 0 senpol 4 -0.0102041 2.4489796 4.4897959 14.4489796 10 0
471 1 1 0 senbre 2 0.9897959 1.4489796 5.4897959 13.4489796 10 0
472 1 1 0 vacnil 7 0.9897959 2.4489796 5.4897959 14.4489796 10 0
473 1 1 1 senbre 2 1.9897959 -2.5510204 6.4897959 9.4489796 10 0
474 1 1 1 albhar 1 1.9897959 -1.5510204 6.4897959 10.4489796 10 0
475 1 1 1 vactor 11 2.9897959 -2.5510204 7.4897959 9.4489796 10 0
476 1 1 1 vacsey 9 2.9897959 -1.5510204 7.4897959 10.4489796 10 0
477 1 1 1 vacnil 7 3.9897959 -2.5510204 8.4897959 9.4489796 10 0
478 1 1 1 senmel 3 3.9897959 -1.5510204 8.4897959 10.4489796 10 0
479 1 1 1 vacsie 10 1.9897959 -0.5510204 6.4897959 11.4489796 10 0
480 1 1 1 senpol 4 1.9897959 0.4489796 6.4897959 12.4489796 10 0
481 1 1 1 vacger 6 2.9897959 -0.5510204 7.4897959 11.4489796 10 0
482 1 1 1 vacrob 8 2.9897959 0.4489796 7.4897959 12.4489796 10 0
483 1 1 1 vacxan 12 3.9897959 -0.5510204 8.4897959 11.4489796 10 0
484 1 1 1 sensen 5 3.9897959 0.4489796 8.4897959 12.4489796 10 0
485 1 0 0 vacxan 12 1.9897959 1.4489796 6.4897959 13.4489796 10 0
486 1 0 0 vactor 11 1.9897959 2.4489796 6.4897959 14.4489796 10 0
487 1 0 0 vacsey 9 2.9897959 1.4489796 7.4897959 13.4489796 10 0
488 1 0 0 vacnil 7 2.9897959 2.4489796 7.4897959 14.4489796 10 0
489 1 0 0 vacger 6 3.9897959 1.4489796 8.4897959 13.4489796 10 0
490 1 0 0 albhar 1 3.9897959 2.4489796 8.4897959 14.4489796 10 0

We can fill in variables to see the model for tree \(i=1\):

\[\begin{align*} y_1 &= \beta_0 \\ &\qquad + \beta^{\text{fire}} * 0 + \beta^{\text{drought}} * 1 + \beta^{\text{herbivory}} * 1 \\ &\qquad + \beta^{\text{fire},\text{drought}} * 0 * 1 + \beta^{\text{fire},\text{herbivory}} * 0 * 1 + \beta^{\text{drought},\text{herbivory}} * 1 * 1 \\ &\qquad + \beta^{\text{species}}_{9} \\ &\qquad + \beta^{\text{species},\text{fire}}_{9} * 0 + \beta^{\text{species},\text{drought}}_{9} * 1 + \beta^{\text{species},\text{herbivory}}_{9} * 1 \\ &\qquad + \beta^{\text{SouthInBlock}} * -4.0 + \beta^{\text{EastInBlock}} *\text{EastInBlock}_i*-2.6 \\ &\qquad + \beta^{\text{South}} * -8.5 + \beta^{\text{East}} *\text{East}_i*-14.6 \\ &\qquad + \beta^{\text{block}}_{1} \\ &\qquad + \epsilon_1 \end{align*}\]

Team activity (9 minutes): What would the model be for tree \(i=1\) if it got no treatment (no fire, no drought, no herbivory) ? Keep the species and plant position variables the same.

An answer

\[\begin{align*} y_1 &= \beta_0 \\ &\qquad + \beta^{\text{species}}_{9} \\ &\qquad + \beta^{\text{SouthInBlock}} * -4.0 + \beta^{\text{EastInBlock}} *\text{EastInBlock}_i*-2.6 \\ &\qquad + \beta^{\text{South}} * -8.5 + \beta^{\text{East}} *\text{East}_i*-14.6 \\ &\qquad + \beta^{\text{block}}_{1} \\ &\qquad + \epsilon_1 \end{align*}\]

We can also write our model for the tree heights in code and use it to simulate fake tree heights:

simulated_data = experimental_design %>% 
  mutate(linear_predictor = beta_0 + 
           beta_fire*Fire +
           beta_drought*Drought +
           beta_herbivory*Herbivory +
           beta_fire_drought*Fire*Drought +
           beta_fire_herbivory*Fire*Herbivory +
           beta_drought_herbivory*Drought*Herbivory +
           beta_species[Species_number] + 
           beta_species_fire[Species_number]*Fire +
           beta_species_drought[Species_number]*Drought +
           beta_species_herbivory[Species_number]*Herbivory +
           beta_SouthInBlock*SouthInBlock +
           beta_EastInBlock*EastInBlock +
           beta_South*South +
           beta_East*East + 
           beta_block[Block],
         y = rnorm(n = n(), mean = linear_predictor, sd = sd_tree))

Simulating from the model, we get a collection of tree heights in centimeters. We can make a histogram of these heights to check that the minimum and maximum are reasonable:

hist(simulated_data$y)

Team activity (8 minutes): What do you notice about these simulated fake tree heights ? What seems reasonable ? What seems unreasonable ?

An answer

Tree heights shouldn’t be negative ! Otherwise things look ok. We could modify the model to prevent negative tree heights. We can do this by transforming the height variable, or by using a model other than the Normal distribution. We will explore this in future courses.

What isn’t in the model ?

Team activity (8 minutes): We could make our mathematical story more complicated. What isn’t in the model ?

An answer

We could add a 3-way interaction between fire, drought, and herbivory. We could let all the interactions vary by species. We could include subplot effects.


More complicated models are more difficult to understand and use. If a simpler model is realistic-enough, we can start there.

Fitting the model

Great ! We have a mathematical story about how the data could be generated:

  • prior: A probability distribution for effects \(p(\beta^{\text{fire}}, \beta^{\text{drought}}, \beta^{\text{herbivory}},...)\). This is our uncertainty about effects before we collect data.
  • data: A probability distribution for tree heights given these effects \(p(y_i\ |\ \beta^{\text{fire}}, \beta^{\text{drought}}, \beta^{\text{herbivory}},...)\).

Now what is our uncertainty about effects after we collect data ? We fit the model (applying Bayes’ Theorem) to get:

  • posterior: A probability distribution for these effects given data on tree heights \(p(\beta^{\text{fire}}, \beta^{\text{drought}}, \beta^{\text{herbivory}},...\ |\ y_i)\)
Bayesian model fitting references We will use an R package called brms, which stands for “Bayesian Regression Models using Stan” (Stan fits Bayesian models).

Fitting to simulated data

Before we have real data, we can learn a lot from fitting the data we simulated from our model.

How much precision do we have to estimate the effect of Fire \(\beta^{\text{fire}}\) ? What about the effect of Drought \(\beta^{\text{drought}}\) ? Herbivory \(\beta^{\text{herbivory}}\) ?

First, recall that our prior for these effects is the same:

\(\beta^{\text{fire}} \sim N(-2,1)\)

\(\beta^{\text{drought}} \sim N(-2,1)\)

\(\beta^{\text{herbivory}} \sim N(-2,1)\)

In other words, prior to data collection, we think fire, drought, and herbivory are both harmful to acacia growth. Our best guess is that they reduce height by 2 centimeters, but we would not be too surprised if it were up to 4 centimeter reduction. We think none help acacia growth.

Then we get information from the data. We combine the prior information with the data information when we fit the model. This gives us posterior guesses about the effects:

calculations to show prior, data, posterior

Let

\(y =\) point in observation space.

\(\tilde{y} =\) observed/realized data, i.e. one such \(y\).

\(a =\) parameter of interest.

\(b =\) other parameter.

\(p(y | a,b) p(a,b) = p(y, a, b) =\) full Bayesian model, with observational model and prior.

To make a prior, likelihood, posterior plot for \(a\), we can get prior \(p(a)\) by the margin of \(p(a,b)\). We can get the posterior \(p(a | \tilde{y})\) by fitting the full Bayesian model and taking the margin of the joint posterior \(p(a, b | \tilde{y})\).

Here is my math to get the marginal likelihood \(p(\tilde{y} | a)\):

Let \(p(y | a,b) p(b | a) * 1 = p_\text{flat}(y, a, b) =\) a modified full Bayesian model with a flat prior on \(a\), but the same \(p(y, b | a)\) as the original full Bayesian model above.

So \(p_\text{flat}(y | a) = p(y | a) = \int p(y, b | a) d b\) is the same as in the original model. So the \(p(\tilde{y} | a)\) that we want can be gotten from \(p_\text{flat}(\tilde{y} | a)\).

Fitting \(p_\text{flat}(y, a, b)\) gives the posterior \(p_\text{flat}(a, b | \tilde{y})\). As usual, we can take the marginal posterior for the parameter of interest to get p_flat(a | ytilde).

By Bayes rule \(p_\text{flat}(a | \tilde{y}) = \frac{p_\text{flat}(\tilde{y} | a) p_\text{flat}(a) }{ p_\text{flat}(\tilde{y})}\) so the posterior \(p_\text{flat}(a | \tilde{y})\) we get is proportional to the marginal likelihood \(p_\text{flat}(\tilde{y} | a)\) that we want.

# keep Fire and Herbivory effect priors separate for an analysis we do below:
priors_no_Fire_or_Herbivory = 
  c(prior(normal(20,5), class = b, coef = Intercept),
    prior(normal(-4,2), class = b, coef = Drought),
    prior(normal(-2,1), class = b, coef = Fire:Drought),
    prior(normal(-2,1), class = b, coef = Fire:Herbivory),
    prior(normal(-2,1), class = b, coef = Drought:Herbivory),
    prior(exponential(0.25), class = sd, coef = Intercept, group = Species),
    prior(exponential(1), class = sd, coef = Fire, group = Species),
    prior(exponential(1), class = sd, coef = Drought, group = Species),
    prior(exponential(1), class = sd, coef = Herbivory, group = Species),
    prior(normal(0,0.05), class = b, coef = SouthInBlock),
    prior(normal(0,0.05), class = b, coef = EastInBlock),
    prior(normal(0,0.05), class = b, coef = South),
    prior(normal(0,0.05), class = b, coef = East),
    prior(student_t(3,0,2.5), class = sigma))

priors = c(priors_no_Fire_or_Herbivory,
           prior(normal(-2,1), class = b, coef = Fire),
           prior(normal(-2,1), class = b, coef = Herbivory))


y_model = brmsformula(y ~ 0 + Intercept +
                        Fire + Drought + Herbivory + 
                        Fire:Drought + Fire:Herbivory + Drought:Herbivory +
                        (1 | Species) + 
                        (-1 + Fire  | Species) + (-1 + Drought | Species) + (-1 + Herbivory | Species) + 
                        SouthInBlock + EastInBlock + 
                        South + East + 
                        (1 | Block))
fit <- brm(y_model,
           prior = priors,
           data = simulated_data
)
priors_Fire_flat = c(priors_no_Fire_or_Herbivory,
                     prior(normal(-2,1), class = b, coef = Herbivory))


priors_Herbivory_flat = c(priors_no_Fire_or_Herbivory,
                          prior(normal(-2,1), class = b, coef = Fire))

fit_Fire_flat <- brm(y_model,
                     prior = priors_Fire_flat,
                     data = simulated_data
)

fit_Herbivory_flat <- brm(y_model,
                          prior = priors_Herbivory_flat,
                          data = simulated_data
)
# Fire
# prior normal(-2,1)
-2 - 2*1
## [1] -4
-2 + 2*1
## [1] 0
beta_fire
## [1] -1.816357
# likelihood
summary(fit_Fire_flat)$fixed["Fire","l-95% CI"]
## [1] -5.284953
summary(fit_Fire_flat)$fixed["Fire","u-95% CI"]
## [1] -0.243377
# posterior
summary(fit)$fixed["Fire","l-95% CI"]
## [1] -3.814237
summary(fit)$fixed["Fire","u-95% CI"]
## [1] -0.8666672
# z-score
(summary(fit)$fixed["Fire","Estimate"] - beta_fire)/summary(fit)$fixed["Fire","Est.Error"]
## [1] -0.7316484
# contraction
1 - summary(fit)$fixed["Fire","Est.Error"]^2/1^2
## [1] 0.4560785
# Herbivory
# prior normal(-2,1)
-2 - 2*1
## [1] -4
-2 + 2*1
## [1] 0
beta_herbivory
## [1] -0.4047192
# likelihood
summary(fit_Herbivory_flat)$fixed["Herbivory","l-95% CI"]
## [1] -1.627991
summary(fit_Herbivory_flat)$fixed["Herbivory","u-95% CI"]
## [1] -0.1244171
# posterior
summary(fit)$fixed["Herbivory","l-95% CI"]
## [1] -1.764698
summary(fit)$fixed["Herbivory","u-95% CI"]
## [1] -0.3497679
# z-score
(summary(fit)$fixed["Herbivory","Estimate"] - beta_herbivory)/summary(fit)$fixed["Herbivory","Est.Error"]
## [1] -1.760362
# contraction
1 - summary(fit)$fixed["Herbivory","Est.Error"]^2/1^2
## [1] 0.8729253
Estimate Est.Error
Intercept 19.8 0.8
Fire -2.4 0.7
Drought -5.9 0.4
Herbivory -1.0 0.4
SouthInBlock 0.0 0.0
EastInBlock 0.0 0.0
South 0.0 0.0
East -0.1 0.0
Fire:Drought -1.6 0.2
Fire:Herbivory -3.1 0.2
Drought:Herbivory -1.3 0.2

Team activity (7 minutes): Which effect are we least certain about in the posterior, \(\beta^{\text{fire}}\), \(\beta^{\text{drought}}\), or \(\beta^{\text{herbivory}}\) ?

An answer

We estimate \(\beta^{\text{fire}}\) to be -1.7 centimeters with estimated error of 0.7 centimeters. We estimate \(\beta^{\text{drought}}\) to be -5.8 centimeters with estimated error of 0.4 centimeters. We estimate \(\beta^{\text{herbivory}}\) to be -0.9 centimeters with estimated error of 0.4 centimeters. So we have least precision for the fire effect.


Improving precision for the fire effect

Team activity (8 minutes): Recall the experimental design. Why do you think precision is worse for the fire effect than for drought and herbivory effects ?

An answer

When we assign fire at the block level, we cannot distinguish between the effect of being in a block versus the effect of fire. For example, if block 3 got fire but block 7 did not get fire, is the difference between the plant heights because of the fire treatment ? Or because block 3 had slightly different soil than block 7 ? Or maybe their rainout shelters are built slightly differently ? When treatment (herbivory and drought) is assigned at a lower level like a subplot, we can compare within the same block and isolate the effect of treatment alone.

But manipulating fire is more difficult than manipulating drought or herbivory. It would be difficult to do the controlled burns at the subplot level.

Team activity (8 minutes): How would you propose to improve the precision of the fire effect ?

An answer

Making the rainout shelters as similar as possible will reduce the size of the block effects. This will make it possible to estimate the effect of fire with more precision. Let’s see this in the mathematical story.

Previously, we assumed block effects with standard deviation 2 centimeters:

\[\beta^{\text{block}}_b \sim N(0,2)\]

Now let’s assume block effects have standard deviation of only 0.2 centimeters:

\[\beta^{\text{block}}_b \sim N(0,0.2)\]

sd_block = 0.2
beta_block = rnorm(n = 10, mean = 0, sd = sd_block)

We refit the model with these new block effects and see these posteriors for the effects:

Estimate Est.Error
Intercept 19.7 0.6
Fire -1.8 0.2
Drought -5.5 0.4
Herbivory -1.3 0.4
SouthInBlock 0.0 0.0
EastInBlock 0.0 0.0
South 0.0 0.0
East -0.1 0.0
Fire:Drought -1.8 0.1
Fire:Herbivory -2.6 0.1
Drought:Herbivory -1.4 0.1
Randomization at block-level vs subplot-level references

Imbens (2014)

10 treatment states scattered about. In this case, assuming iid errors does not seem right: instead, one might cluster the errors by state. Now when above a certain amount increasing within-state sample size won’t do much for power…

Suppose we think of the outcomes being affected by individual level characteristics as well as state effects, with the latter modeled as random effects, independent of the treatment…the variance estimates based on assuming there are no state effects would be underestimating the true variance

Gelman et al. (2022)

As Imbens (2014) points out, with only one treatment village and one control village, the treatment e!ect cannot be separated from the difference in village effects.

Su and Ding (2021)

Cluster-randomized experiments are widely used due to their logistical convenience and policy relevance. To analyse them properly, we must address the fact that the treatment is assigned at the cluster level instead of the individual level…Feller and Gelman (2015) reviewed the standard approach of using hierarchical models for causal inference with cluster-randomized experiments.

Feller and Gelman (2015)

In cluster-randomized experiments, every unit in the cluster receives the same treatment; in other words, randomization occurs at the group level…

Fitting to real data

We will have a follow-up course to make this step !

Team activity (7 minutes): Most experiments are not implemented exactly how they are designed. We need to model what was actually implemented, including any logistical problems or changing strategies. What could have changed from the design with this experiment ? Ask Arjun and Charles !

An answer

Cohort effects: We did 2 rounds of seed planting. We germinated seeds at start and again at the end of rainy season. So for all species there are two age groups. This increases the variation in size.

Subplot layout: To fit 6 plant positions within a square subplot, we used this layout from wikipedia:

We have the true layout coordinates and will use those in the analysis.

The corn plant subplots were implemented instead as low-density subplots, with corn in plant position 6, vacxan acacia in position 3, and a randomly-chosen acacia species in position 2. The species was sampled without replacement, so 10 of the 12 species are represented as the control plant (one for each block).


Baseline data

We have baseline (pre-treatment) data on plant heights. We can use this data to look at differences across treatment groups, as suggested by Imai et al. (2008):

Researchers analysing data from randomized experiments that did not block on all observed covariates can check balance, but they do not need hypothesis tests to do so….In any study where all observed covariates were not fully blocked ahead of time, balance should be checked routinely by comparing observed covariate differences between the treated and control groups

In the analysis of the endline (post-treatment) data, we can use the baseline (pre-treatment) data to improve our estimates.

In ecology this is sometimes called Before-After-Control-Impact (BACI), see for example Chevalier et al. (2019).

Statisticians advocate for this, for example Gelman et al. (2019):

The treatment and placebo groups differed in their pre-treatment levels of exercise time, with mean values of 528.0 and 490.0 s, respectively. This sort of difference is fine - randomization assures balance only in expectation - but it is important to adjust for this discrepancy in estimating the treatment effect.

Celebrate

Designing a study is hard and we deserve to celebrate before we even collect data.

Team activity (23 minutes): What went right in this experimental design ? What do you like about it ? How will you celebrate these design ideas in your own work ?

An answer

My favorite aspects of this experimental design:

  • treatment implementation and outcome measurements are designed to be close to the real-world aspects we want to study
  • control groups allow for estimating effects of treatments
  • randomization allows for estimating these effects without bias
  • careful setup to make groups comparable, e.g. building rainout shelters as similarly as possible, which helps improve precision
  • data collection from baseline until the endline of the study enables many analyses, a full time series ! Helps us learn a lot and be more efficient with the data by adjusting for baseline differences, see e.g. Gelman et al. (2019).