-
Notifications
You must be signed in to change notification settings - Fork 6
Expand file tree
/
Copy pathgxpopbehaviorhabitatq1.Rmd
More file actions
4814 lines (3743 loc) · 335 KB
/
gxpopbehaviorhabitatq1.Rmd
File metadata and controls
4814 lines (3743 loc) · 335 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: "Implementing a rapid geographic range expansion - the role of behavior changes"
author:
- '[Logan CJ](http://CorinaLogan.com)^1^*'
- '[McCune KB](https://www.kelseymccune.com/)^2^'
- "LeGrande-Rolls C^1^"
- Marfori Z^1^
- Hubbard J^3^
- '[Lukas D](http://dieterlukas.mystrikingly.com)^1^'
output:
html_document:
toc: yes
toc_depth: 5
toc_float:
collapsed: no
code_folding: hide
md_document:
toc: yes
word_document:
toc: yes
toc_depth: '5'
pdf_document:
keep_tex: yes
latex_engine: xelatex
bibliography: MyLibrary.bib
csl: Peer-Community-Journal-PCI.csl
urlcolor: blue
header-includes:
- \usepackage{pdflscape}
- \newcommand{\blandscape}{\begin{landscape}}
- \newcommand{\elandscape}{\end{landscape}}
- \usepackage{longtable} #set longtable = TRUE in kable() to have a long table break over multiple pages https://bookdown.org/yihui/bookdown/tables.html
#For posting final version at EcoEvoRxiv: replace pre-introduction with https://docs.google.com/document/d/1DtcfJ9zE6c5MawGJ2Uq94Pmav-TldhTTKDmaGMJYJ7I/edit?usp=sharing . The PCI templates are located at https://osf.io/hfdjc
---
```{r setup, include=FALSE}
library(knitr)
library(formatR)
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=70),tidy=TRUE)
#Make code chunks wrap text so it doesn't go off the page when knitting to PDF
knitr::opts_chunk$set(echo=F, include=T, results='asis', warning=F, message=F)
#sets global options to display code along with the results https://exeter-data-analytics.github.io/LitProg/r-markdown.html
#set echo=F for knitting to PDF (hide code), and echo=T for knitting to HTML (show code)
```
```{r cleanbib, include=FALSE, eval=F}
### make a bibtex file that has only the references cited in this rmd
#load the Rmd file
Rmd <- readChar("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/gxpopbehaviorhabitatq1.Rmd",nchars=1e7)
#find all in text citations that start with @, but are preceded by a space
pattern <- "\\ @(.*?)\\ "
m <- regmatches(Rmd,gregexpr(pattern,Rmd))[[1]]
m
res<- gsub("\\ ","",m) #delete spaces
res<- gsub("\\]","",res) #delete ]
res<- gsub("\\;","",res) #delete ;
res<- gsub("\\,","",res) #delete ,
res<- gsub("\\.","",res) # delete .
#find all in text citations that start with @, but are preceded by a "["
pattern2 <- "\\[@(.*?)\\ "
m2 <- regmatches(Rmd,gregexpr(pattern2,Rmd))[[1]]
m2
res2<- gsub("\\[","",m2)
res2<- gsub("\\]","",res2)
res2<- gsub("\\]","",res2)
res2<- gsub("\\;","",res2)
res2<- gsub("\\,","",res2)
res2<- gsub("\\.","",res2)
res2<- gsub("\\ ","",res2)
#combine both patterns
allbibtexkeys<-c(res,res2)
#write to a new file and then clean it up manually
write(allbibtexkeys,file="gxpopbehaviorhabitatq1_bibtexkeys.txt")
#load the cleaned txt file
allbibtexkeys<-read.csv("gxpopbehaviorhabitatq1_bibtexkeys.txt")
#use bib2df to convert the bibliography file into a dataframe
install.packages("bib2df")
library(bib2df)
#load the bib from GitHub
df <- bib2df("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/MyLibrary.bib")
#remove the @ to match the entry in the bib file
allbibtexkeys2<- gsub("\\@","",allbibtexkeys[,1])
#filter the full bib to only keep the entries that are cited here
df_filtered<-df[df$BIBTEXKEY %in% allbibtexkeys2,]
#use bib2df to convert the data frame into the bibliography file that only contains the citations in this Rmd
df2bib(df_filtered, file = "xpopq1refs.bib", append = FALSE)
```
[](https://doi.org/10.24072/pci.ecology.100535)
**Affiliations:** 1) Max Planck Institute for Evolutionary Anthropology, 2) University of California Santa Barbara, 3) Animal Behavior Graduate Group, University of California Davis. *Corresponding author: corina_logan@eva.mpg.de
**This research article has been peer reviewed and recommended by:**
Esther Sebastián González (2023) Behavioral changes in the rapid geographic expansion of the great-tailed grackle. *Peer Community in Ecology*, 100535. [https://doi.org/10.24072/pci.ecology.100535](https://doi.org/10.24072/pci.ecology.100535). Reviewers: Francois-Xavier Dechaume-Moncharmont, Pizza Ka Yee Chow, and 1 anonymous reviewer
**Cite as:**
Logan CJ, McCune KB, LeGrande-Rolls C, Marfori Z, Hubbard J, Lukas D. 2023.Implementing a rapid geographic range expansion - the role of behavior changes. EcoEvoRxiv, version 3, peer reviewed and recommended by *Peer Community in Ecology*. doi: [https://doi.org/10.32942/X2N30J](https://doi.org/10.32942/X2N30J)
**This article began as a preregistration, which was pre-study peer reviewed and received an In Principle Recommendation of the [version](https://github.com/corinalogan/grackles/blob/master/Files/Preregistrations/gxpopbehaviorhabitatPassedPreStudyPeerReview16Dec2021.pdf) on 16 Dec 2021 by:**
Esther Sebastián González (2020) The role of behavior and habitat availability on species geographic expansion. *Peer Community in Ecology*, 100062. [10.24072/pci.ecology.100062](https://doi.org/10.24072/pci.ecology.100062). Reviewers: Caroline Nieberding, Tim Parker, and Pizza Ka Yee Chow
## ABSTRACT
It is generally thought that behavioral flexibility, the ability to change behavior when circumstances change, plays an important role in the ability of species to rapidly expand their geographic range. Great-tailed grackles (*Quiscalus mexicanus*) are a social, polygamous species that is rapidly expanding its geographic range by settling in new areas and habitats. They are behaviorally flexible and highly associated with human modified environments, eating a variety of human foods in addition to foraging on insects and on the ground for other natural food items. They offer an opportunity to assess the role of behavior change over the course of their expansion. We compared behavior in wild-caught grackles from two populations across their range (an older population in the middle of the northern expansion front: Tempe, Arizona, and a more recent population on the northern edge of the expansion front: Woodland, California) to investigate whether certain behaviors (flexibility, innovativeness, exploration, and persistence) have higher averages and variances in the newer or older population. We found that grackles in the edge population had a higher flexibility variance (measured by reversal learning) and a higher persistence average (they participated in a larger proportion of trials), and that there were no population differences in average levels of flexibility, innovativeness (number of loci solved on a multiaccess box), or exploration (latency to approach a novel environment). Our results elucidated that individuals differentially expressing a particular behavior in an edge population could facilitate the rapid geographic range expansion of great-tailed grackles, and we found no support for the importance of several traits that were hypothesized to be involved in such an expansion. Our findings highlight the value of population studies and of breaking down cognitive concepts into direct measures of individual abilities to better understand how species might adapt to novel circumstances.
\newpage
## INTRODUCTION
It is generally thought that behavioral flexibility (hereafter, “flexibility”) plays an important role in the ability of a species to rapidly expand their geographic range [e.g., @lefebvre1997feeding; @griffin2014innovation; @chow2016practice; @sol2000behavioural; @sol2002behavioural; @sol2005big; @sol2007big]. It is predicted that flexibility, the ability to change behavior when circumstances change through packaging information and making it available to other cognitive processes [see @mikhalevich_is_2017 for theoretical background on our flexibility definition], as well as exploration (latency to explore a novel environment or object) and innovation [creating new behaviors or using existing behaviors in a new context, @griffin2014innovation] facilitate the expansion of individuals into completely new areas. However, the role of these behaviors in the process of establishing a population in a particular area is predicted to diminish after a certain number of generations [@wright2010behavioral]. In support of this, experimental studies have shown that latent abilities are primarily expressed in a time of need [e.g., @taylor2007spontaneous; @bird2009insightful; @manrique2011spontaneous; @auersperg2012spontaneous; @laumer2018spontaneous].
To determine whether a behavior (e.g., flexibility, innovativeness, exploration, persistence) is involved in a rapid geographic range expansion, direct measures of behaviors in individuals must be collected in populations across the range of the species (see the discussion on the danger of proxies of flexibility in Logan et al., 2018). Flexibility, the ability to recognize that something about the environment has changed and decide to consider other options for deploying behavior [@mikhalevich_is_2017], is distinct from innovation, which is the specific stringing together of particular actions in a new way or in a new context [@griffin2014innovation]. Innovative behavior can be related to flexibility in that it may occur in response to the decision to change behavior in some way. Investigations of behavior in invasive species and species that are rapidly expanding their geographic ranges that compare edge versus core populations are rare. Behavioral evidence from invasive species indicates that Common mynas (*Sturnus tristis*) on the invasion front are more innovative than those from populations away from the front as well as those in their native range [@cohen2020innovation]. Similarly, spiders (*Cyrtophora citricola*) and bank voles (*Myodes glareolus*) from edge populations are less exploratory than those from core populations [@chuang2021personality; @eccard2022timid]. An increase in innovation in newly established populations could facilitate new foraging techniques and the ability to exploit new food sources [@griffin2016invading], while a decrease in exploration could reduce their risk of encountering danger in a new area. More data from more species is required to discover whether these results are generalizable to an invasion or rapid range expansion context. As such, it is important to decide which measures are the best proxies of the behavior in question. For example, exploration is often measured as activity levels [e.g., @fox_behavioural_2009; @logan2016behavioral], however it is important to distinguish activity levels, which could be an indicator of stress, from the curiosity to investigate novelty [@mettke2002significance]. The latter can be accomplished by placing a novel environment or object inside of the familiar environment, thus making it optional to approach the novel element. Additionally, we can distinguish exploration from boldness through variation in food deprivation or placement of food. For boldness, the behavioral response to a potential threat, subjects are usually food deprived and then a preferred food item is placed next to the novel object [@reale2007integrating]. Whereas, in exploration assays, the regular maintenance diet is provided far away from the novel element to assess the willingness to investigate novelty without the need for food [@mettke2002significance]. The latter ensures that the individual approaches the novel element primarily because they are internally motivated to explore something new.
Persistence behavior, “a measure of task-directed motivation” (Griffin & Guez 2014), to our knowledge, has not been investigated across populations of species that are rapidly expanding their geographic ranges. However, it could facilitate a range expansion through improving problem solving success [@morand2011innovators] and efficiency [@chow2016practice]. There is some indication that this could be the case in a cross-species comparison of Invasive mynas who were found to be more persistent than native noisy miners (*Manorina melanocephala*) even though both species are successful in urban environments [@griffin2015innovative]. Persistence is measured in a variety of ways [e.g., work time, number of touches to the test apparatus, number/frequency of unsuccessful manipulations, etc., see @griffin2014innovation for a review], which makes it a difficult variable to compare across studies. Many measures of persistence are resource intensive to collect because they involve hundreds of hours of video coding, which could prohibit some researchers from being able to measure this variable due to time and financial constraints. Therefore, we developed an easy to calculate measure that we believe better represents task-directed motivation in grackles: the number of trials participated in divided by the total number of trials offered.
We expect that the actual act of continuing a range expansion relies on flexibility, exploration, innovation, and persistence. It is therefore likely that these behaviors are expressed more on the edge of the expansion range where there have not been many generations to accumulate relevant knowledge about or genetic adaptations to the environment. Our study aims to test whether behavioral flexibility, innovativeness, exploration, and persistence play a role in the rapid geographic range expansion of great-tailed grackles (*Quiscalus mexicanus*). Great-tailed grackles are behaviorally flexible [@logan2016flexibilityproblem], rapidly expanding their geographic range [@wehtje2003range], and highly associated with human modified environments [@johnson2001great], thus offering an opportunity to assess the role of behavior across their expansion. This social, polygamous species eats a variety of human foods in addition to foraging on insects and on the ground for other natural food items [@johnson2001great]. This opportunistic foraging behavior increases the ecological relevance of comparative cognition experiments that measure individual behavior abilities: grackles eat at outdoor cafes, from garbage cans, and on crops, where they generally gain experience in the wild with approaching and opening novel objects to seek food (e.g., attempting to open a ketchup packet at an outdoor cafe, climbing into garbage cans to get french fries at the zoo, dunking sugar packets in water). Consequently, tests involving human-made apparatuses are ecologically relevant for this species. We compared behavior in wild-caught great-tailed grackles from two populations across their range. We use previously published data from @logan2023flexmanip for an older population in the middle of the northern expansion front in Tempe, Arizona, as well as new data collected on a more recent population on the northern edge of the expansion front in Woodland, California (Figure 1, Table 1). We investigated whether certain behaviors had higher averages and variances in the edge population relative to the older population. Specifically, we investigated behavioral flexibility, measured as reversal learning of food-filled colored tube preferences [@logan2016behavioral; @logan2023flexmanip]; innovativeness, measured as the number of loci they solve to access food from a puzzle box [@auersperg_flexibility_2011; @logan2023flexmanip]; exploration, measured as the latency in seconds to approach a novel environment in the absence of nearby food [@mccune2019exploration; @mettke2009spatial]; and persistence, measured as the proportion of trials they participated in during the flexibility and innovativeness experiments (Figure 2). While it is possible for individuals in the wild to learn asocially and socially about new foods or foraging techniques to assess whether the risks are low enough to encourage exploration behavior, we focused on measuring these four behaviors in an asocial context to allow us to obtain the individual’s actual cognitive performance (i.e., in the absence of dominant individuals who might hinder subordinates from participating). There could be multiple mechanisms underpinning the results, however our aim was to narrow down the role of changes in behavior in the range expansion of great-tailed grackles.

\newpage
**Table 1.** Population characteristics for the field sites. The number of generations at a site is based on a generation length of 5.6 years for this species [@GTGRbirdlife2018, note that this species starts breeding at age 1] and on the first year in which this species was reported (or estimated) to breed at each location (Woodland, California: Yolo Audubon Society’s newsletter *The Burrowing Owl* from July 2004; and Tempe, Arizona: estimated based on 1945 first-sighting report in nearby Phoenix, Arizona [@wehtje2004thesis] to which we added 6 years to account for the average time between first-sighting and first-breeding - see Table 3 in @wehtje2003range. The average number of generations was calculated using the number of years of breeding (the “Breeding since” year up to 2020, the final year of data collection in Tempe, and 2022, the final year of data collection in Woodland) divided by the 5.6 year generation length.
```{r table1, eval=TRUE, warning=FALSE, results='asis', include=TRUE}
d <- read.csv(url("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/gxpopbehaviorhabitatq1_Table1.csv"), header=F, sep=",", stringsAsFactors=F)
colnames(d) <- c("Site","Range position","Breeding since","Number of years breeding","Average number of generations","Citation")
library(kableExtra)
knitr::kable(d) %>%
kable_styling(full_width = T, position = "left")
```
### RESEARCH QUESTION: Are there differences in behavioral traits (flexibility, innovation, exploration, and persistence) between populations across the great-tailed grackle's geographic range?
**Prediction 1:** **If behavior modifications are needed to adapt to new locations, then there is a higher average and/or larger variance of at least some traits (behaviors) thought to be involved in range expansions** (behavioral flexibility: speed at reversing a previously learned color preference based on it being associated with a food reward; innovativeness: number of options solved on a puzzle box; exploration: latency to approach/touch a novel object; and persistence: proportion of trials participated in with higher numbers indicating a more persistent individual) **in the grackles sampled from the more recently established population relative to the individuals sampled in the older population** (Table 1). Higher **averages** in behavioral traits indicate that each individual can exhibit more of that trait (e.g., they are more flexible/innovative/exploratory/persistent). Perhaps in newly established populations, individuals need to learn about and innovate new foraging techniques or find new food sources. Perhaps grackles require flexibility to visit these resources according to their temporal availability and the individual's food preferences. Perhaps solving such problems requires more exploration and persistence. Higher **variances** in behavioral traits will indicate that there is a larger diversity of individuals in the population, which means that there is a higher chance that at least some individuals in the population could innovate foraging techniques and be more flexible, exploratory, and persistent, which could be learned by conspecifics and/or future generations. *This supports the hypothesis* that changes in behavioral traits facilitate the great-tailed grackle's geographic range expansion.
## METHODS
### Sample
Great-tailed grackles were caught in the wild in Woodland and in the Bufferlands of Sacramento, California. Some of our banded individuals were found in Woodland and the Bufferlands, which are 32 km apart, therefore we considered this one population. We caught grackles with walk-in traps and mist nets. Mist nets decrease the likelihood of a selection bias for exploratory and bold individuals because grackles cannot see the trap. We aimed to bring adult grackles, rather than juveniles, temporarily into the aviaries for behavioral choice tests to avoid the potential confound of variation in cognitive development due to age, as well as potential variation in fine motor-skill development [e.g., early-life experience plays a role in the development of holding/grasping objects, @collias1964development; @rutz2016discovery] with variation in our target variables of interest. Observations from members of the Yolo Audubon Society in Woodland, Davis, and Sacramento, California suggest that movement into new areas is most likely by adults or groups of mixed age individuals (Yolo Audubon Society’s newsletter *The Burrowing Owl*). Accordingly, if there are differences associated with presence at the edge of the rane, these differences should also be expressed in adults. Adults were identified from their eye color, which changes from brown to yellow upon reaching adulthood [@johnson2001great]. However, due to difficulties in trapping this species at this site, we also tested some juveniles. This should not pose a problem because we found that the two juveniles (Taco and Chilaquile) we tested in the Tempe population did not perform differently from adults [@logan2023flexmanip; @blaisdell2021causal; @logan2021inhibition; @seitz2021touchscreentraining]. We applied colored leg bands in unique combinations for individual identification. Some individuals (n=33 in Woodland) were brought temporarily into aviaries for behavioral choice tests, and then released back to the wild at their point of capture. Grackles were individually housed in an aviary (each 244 cm long by 122 cm wide by 213 cm tall) for 3 weeks to 6 months where they had *ad lib* access to water at all times and were fed Mazuri Small Bird maintenance diet *ad lib* during non-testing hours (minimum 20 h per day), and various other food items (e.g., peanuts, bread) during testing (up to 4 h per day per bird). Individuals were given three to four days to habituate to the aviaries and then their test battery began on the fourth or fifth day (birds were usually tested six days per week, therefore if their fourth day occurred on a day off, they were tested on the fifth day instead).
We tested as many great-tailed grackles as we could during the 2 years we spent at each of our field sites given that the birds were only brought into the aviaries during the non-breeding season (September through April). It is time intensive to conduct the aviary test battery (3 weeks-6 months per bird), therefore we aimed to meet the minimum sample sizes in Supplementary Material 1 and 2. We aimed for an equal sex ratio of subjects (50% female) and achieved an overall 47% female (this percentage differed depending on the test). We expected to test 20 grackles per site. See the gxpopbehaviorhabitat_data_testhistory.csv data sheet at @logan2023xpopdata for a list of the order of experiments for each individual at the Woodland site, and g_flexmanip_data_AllGrackleExpOrder.csv at @logan2023flexmanipdata for the Tempe grackles.We stopped collecting data on wild-caught great-tailed grackles once we met our minimum sample size (Supplementary Material 1 and 2).
### Protocols
Experimental and habituation protocols are available in Supplementary Material 5. In brief, the **flexibility** protocol [from @logan2023flexmanip] used reversal learning with color tubes. Grackles were first habituated to a yellow tube and trained to search for hidden food. A light gray tube and a dark gray tube were placed on the table or floor: one color always contained a food reward (not visible by the bird) while the other color never contained a reward. The bird was allowed to choose one tube per trial. An individual was considered to have a preference if it chose the rewarded option at least 85% of the time (17/20 correct) in the most recent 20 trials (with a minimum of 8 or 9 correct choices out of 10 on the two most recent sets of 10 trials). We used a sliding window in 1-trial increments to calculate whether they passed after their first 20 trials. Once a bird learned to prefer one color, the contingency was reversed: food was always in the other color and never in the previously rewarded color. The flexibility measure was how many trials it took to reverse their color preference using the same passing criterion. The first rewarded color in reversal learning was counterbalanced across birds. The rewarded option was pseudorandomized for side (and the option on the left was always placed on the substrate first by the experimenter). Pseudorandomization consisted of alternating location for the first two trials of a session and then keeping the same color on the same side for at most two consecutive trials thereafter. A list of all 88 unique trial sequences for a 10-trial session, following the pseudorandomization rules, was generated in advance for experimenters to use during testing (e.g., a randomized trial sequence might look like: LRLLRRLRLR, where L and R refer to the location, left or right, of the rewarded tube). Randomized trial sequences were assigned randomly to any given 10-trial session using a random number generator (random.org) to generate a number from 1-88.
The **innovativeness** protocol [from @logan2023flexmanip; and based on the experimental design by @auersperg_flexibility_2011] used a multiaccess log. Grackles were first habituated to the log apparatus with all of the doors locked open and food inside each locus. After habituation, the log, which had four ways of accessing food (pull drawer, push door, lift door up, swing door out), was placed on the ground and grackles were allowed to attempt to solve or successfully solve one option per trial. Once a bird successfully solved an option three times, it became non-functional (the door was locked open and there was no food at that locus). The experiment ended when all four loci became non-functional, if a bird did not come to the ground within 10 min in three consecutive test sessions, or if a bird did not obtain the food within 10 min (or 15 min if the bird was on the ground at 10 min) in three consecutive test sessions.
**Persistence** was measured as the proportion of trials participated in during the flexibility and innovativeness experiments (after habituation, thus it is not confounded with boldness). The higher the number, the more persistent they were. This measure indicates that those birds who do not participate as often were less persistent in engaging with the task. We generally offered a grackle the chance to participate in a trial for 5 min. If they did not participate within that time, we recorded -1 in the data sheet, the apparatus was removed, and the trial was re-attempted later.
**Exploration** was measured as the latency to approach within 20 cm of a novel environment inside of their familiar aviary environment and this test was conducted two times for each bird so we could obtain individual consistency measures. Time 1 occurred on the individual's 8th day in the aviary and Time 2 occurred 1 week after Time 1. The bird’s regular food was moved to one end of the aviary, away from the novel environment, and we first conducted a motivation test where we placed a piece of preferred food on the ground and waited out of view for 5 min. We only proceeded with the exploration assay if the bird ate the food. This motivation test allowed us to determine whether the grackle was interested in coming to the ground at all, where, for example, a grackle might not eat the food because it has just bathed and is primarily focused on preening and drying feathers. The bird was then exposed to first a familiar environment without the novel environment for 45 min and then to a novel environment (a tent) that is placed on the ground within the familiar environment for 45 min. If an individual did not approach within 20 cm, it was given a latency of 2701 sec (45 min plus 1 sec). In a previous experiment [@mccune2019exploration], we validated that grackles did not perceive the novel environment as threatening (i.e., it was not a measure of boldness).
**Experimental order:** The order of experiments for reversal learning or multiaccess log, was counterbalanced across birds for the Woodland population. The Arizona population received the reversal learning experiment first because their flexibility was manipulated to determine whether manipulating flexibility influences performance on subsequent tests [see @logan2023flexmanip].
![Experimental protocol (see Supplementary Material 5 for more details). Great-tailed grackles from the older and newer populations were tested for their: (top left) flexibility as the number of trials to reverse a previously learned color tube-food association; (middle) innovativeness as the number of loci (lift, swing, pull, push) solved to obtain food from within a multiaccess log; (bottom left) persistence as the proportion of trials participated in during flexibility and innovativeness tests; and (far right) exploration as the latency to approach a novel environment placed inside of the familiar environment with regular food present, but not near the novel environment. The order of the flexibility and innovativeness experiments was counterbalanced for the California grackles and they received their first exploration assay as close as possible to day 8 in the aviaries. The Arizona grackles received the flexibility experiment first (because they underwent a flexibility manipulation) and the innovativeness experiment and exploration assay afterward [note that there could have been other experiments between the flexibility experiment and the innovation experiment and exploration assay because their test battery was much larger than that of the California birds, @logan2023flexmanip]. See the test history for each bird in the gxpopbehaviorhabitatq1_data_testhistory.csv data sheet at @logan2023xpopdata.](gxpopbehaviorhabitatFig2.png)
### Statistical analyses
We used **simulations** and designed **customized models** to determine what sample sizes would allow us to detect differences between sites following the methods in @rethinking2020 [Supplementary Material 1 and 2; see chapter 5.3 in @bolker2008ecological for why simulations perform more powerful power analyses]. We did not **exclude** any data, and if data were **missing** (e.g., if a bird participated in only one of the two experiments) for an individual in a given experiment, then this individual was not included in that analysis. Analyses were conducted in R [current version `r getRversion(), @rcoreteam] and Stan [version 2.18, @carpenter2017stan] using the following packages: psych [@psych] and irr [@gamer2012package] for calculating interobserver reliability scores; rethinking [@rethinking2020], cmdstanr [@cmdstanr], rstan [@rstan], posterior [@posterior] and Rcpp [@rcpp] for conducting Bayesian analyses; knitr [@xie2018knitr; @xie2017dynamic; @xie2013knitr], formatR [@formatr], dplyr [@dplyr], tidyr [@tidyr], kableExtra [@kableextra], lattice [@Sarkar2008], and gridExtra [@gridExtra] for formatting; DHARMa [@hartig2019dharma] for data cleaning; lme4 [@lme4; @bates2012lme4] and MCMCglmm [@hadfield2010mcmc] for running GLMMs; and rptR [@stoffel2017rptr] for calculating repeatability. Interobserver reliability scores indicated high agreement across coders for all dependent variables (see Supplementary Material 3 for details).
##### Flexibility analyses
\
*Model and simulation*
We modified the reversal learning Bayesian model in @blaisdell2021causal to simulate and analyze population differences in reversal learning, and calculated our ability to detect differences between populations. The model accounts for every choice made in the reversal learning experiment and updates the probability of choosing either option after the choice is made depending on whether that choice contains a food reward or not. It does this by updating three main components for each choice: an attraction score (how much an individual prefers one option over the other), a learning rate ($\phi$; higher values mean the individual updates their attraction score at a higher rate), and a rate of deviating from learned attractions ($\lambda$; lower values mean the individual is choosing between the options more randomly). The attraction score is the weight an individual gives to a particular option based on its past reward history for that option with attractions increasing if they received a reward when previously choosing that option. The decision regarding which of the two options to make is determined by the relative weights of the two attraction scores (each option gets its own attraction score).
As in @blaisdell2021causal, we, too, used previously published data on reversal learning of color tube preferences in great-tailed grackles in Santa Barbara, California [@logan2016behavioral] to inform the model modifications. We modified the @blaisdell2021causal model in a two ways: 1) we set the initial attraction score assigned to option 1 and option 2 (the light gray and dark gray tubes) to 0.1 rather than 0.0 [see @lukasflexmanip for more detail]. This change assumes that there would be some inclination (rather than no inclination) for the bird to approach the tubes when they are first presented because they are previously trained to expect food in tubes. This also allows the attraction score to decrease when a non-rewarded choice is made near the beginning of the experiment. With the previous initial attraction scores set to zero, a bird would be expected to choose the rewarded option in 100% of the trials after the first time it chose that option (attraction cannot be lower than zero, and choice is shaped by the ratio of the two attractions so that when one option is zero and the other is larger than zero, the ratio will be 100% for the rewarded option). 2) We changed the updating so that an individual only changes the attraction toward the option they chose in that trial [either decreasing their attraction toward the unrewarded option or increasing their attraction toward the rewarded option; see @lukasflexmanip for more detail]. Previously, both attractions were updated after every trial, assuming that individuals understand that the experiment is set up such that one option is always rewarded. For our birds, we instead assumed that individuals will focus on their direct experience rather than making abstract assumptions about the test. Our modification resulted in needing a higher $\phi$ to have the same learning rate as a model where both attraction scores update after every trial [@lukasflexmanip]. This change also appears to better reflect the performance of the Santa Barbara grackles, because they had higher $\phi$ values, which, in turn, meant lower $\lambda$ values to reflect the performance during their initial learning. These lower $\lambda$ values better reflected the birds' behavior during the first reversal trials: a large $\lambda$ value means that birds continue to choose the now unrewarded option almost 100% of the time, whereas the lower $\lambda$ values mean that birds start to explore the rewarded option relatively soon after the switch of the rewarded option [@lukasflexmanip].
We first reanalyzed the Santa Barbara grackle data to obtain the $\phi$ and $\lambda$ values with this revised model, which informed our expectations of what a site’s mean and variance might be. Then we ran simulations, where we determined that we wanted to make the previously mentioned modifications to the stan model. This model was used to analyze the actual data after it was collected, using only data from the first reversals to eliminate the need to modify the model to include treatment (whether an Arizona grackle was manipulated or not). We used an analysis called a contrast to assess whether one site was systematically larger or smaller than the other by estimating what percentage of each sample of differences is either larger or smaller than zero. If 89% of the differences are larger than zero, then the older population has a larger mean, and if 89% of the differences are smaller than zero, then the edge population has a larger mean. If 89% of the differences cross zero, then we conclude that there is no strong difference between the sites. See Supplementary Material 1 and 2 for more details. To determine whether there were differences between the variances in $\phi$ and $\lambda$ between sites, we conducted models as follows:
$\phi_{i}$ or $\lambda_{i}$ ~ Normal($\mu$, $\sigma$[site]) [likelihood],
log($\mu$) ~ $\alpha$[site] [model],
where either $\phi_{i}$ or $\lambda_{i}$ were used as the response variable, $\sigma$[site] allows a separate variance to be assigned to each site, $\alpha$ is the intercept for the $\phi_{i}$ or $\lambda_{i}$ means, and each site gets its own intercept. We then ran a contrast to determine whether there was a difference in variances between the sites.
```{r sim_power_philambda, eval=F, warning=FALSE, results='asis', include=TRUE}
library(rethinking)
library(dplyr)
# We first reanalyze the Santa Barbara data to obtain values for phi and lambda upon which to base our expectations of what the mean and variance might be at a site.
# We load the Santa Barbara data from github
d <- read.csv(url("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/g_causal_data_SBreverse.csv"), header=T, sep=",", stringsAsFactors=F)
# We need to perform some modifications on these data to set them up for the STAN model
d$id <- sapply(1:nrow(d), function (i) which(unique(d$Bird) == d$Bird[i]) )
d <- d[-which(d$Experiment=="Refresher"),]
d <- d[with(d, order(d$id)), ]
d$Choice <- NA
for (i in 1: nrow(d)) {
if (d$Experiment[i] == "Initial"){
if (d$Correct[i] == 1){
d$Choice[i] <- 1
} else {
d$Choice[i] <- 2
}
} else {
if (d$Correct[i] == 1){
d$Choice[i] <- 2
} else {
d$Choice[i] <- 1
}
}
}
d[d=="Initial"] <- 0
d[d=="Reverse"] <- 1
d$Experiment <- as.integer(d$Experiment)
# We can now extract the relevant data for the STAN model. In the model, we want to estimate the phi and lambda for each individual that was tested in Santa Barbara based on their choices during the initial learning phase and the reversal learning phase.
dat <- as.list(d)
dat$N <- nrow(d)
dat$N_id <- length(unique(d$id))
# The STAN model is set up to have the inital attraction for each option set to 0.1, and that individuals only learn the reward of the option they chose in a given trial.
reinforcement_model_nonzeroattraction <- "
data{
int N;
int N_id;
int id[N];
int Trial[N];
int Choice[N];
int Correct[N];
}
parameters{
real logit_phi;
real log_L;
// Varying effects clustered on individual
matrix[2,N_id] z_ID;
vector<lower=0>[2] sigma_ID; //SD of parameters among individuals
cholesky_factor_corr[2] Rho_ID;
}
transformed parameters{
matrix[N_id,2] v_ID; // varying effects on stuff
v_ID = ( diag_pre_multiply( sigma_ID , Rho_ID ) * z_ID )';
}
model{
matrix[N_id,2] A; // attraction matrix
logit_phi ~ normal(0,1);
log_L ~ normal(0,1);
// varying effects
to_vector(z_ID) ~ normal(0,1);
sigma_ID ~ exponential(1);
Rho_ID ~ lkj_corr_cholesky(4);
// initialize attraction scores
for ( i in 1:N_id ) {
A[i,1] = 0.1; A[i,2] = 0.1';
}
// loop over Choices
for ( i in 1:N ) {
vector[2] pay;
vector[2] p;
real L;
real phi;
// first, what is log-prob of observed choice
L = exp(log_L + v_ID[id[i],1]);
p = softmax(L*A[id[i],1:2]' );
Choice[i] ~ categorical( p );
// second, update attractions conditional on observed choice
phi = inv_logit(logit_phi + v_ID[id[i],2]);
pay[1:2] = rep_vector(0,2);
pay[ Choice[i] ] = Correct[i];
A[ id[i] , Choice[i] ] = ( (1-phi)*(A[ id[i] , Choice[i] ]) + phi*pay[Choice[i]])';
}//i
}
"
# We run this model for the Santa Barbara data
m_SB <- stan( model_code = reinforcement_model_nonzeroattraction, data=dat ,iter = 5000, cores = 4, chains=4, control = list(adapt_delta=0.9, max_treedepth = 12))
# We extract the estimated phi and lambda values for each bird from the STAN model
s <- extract.samples(m_SB)
santabarbara_lambda <- sapply(1 : dat$N_id, function(x) exp( mean(s$log_L) + mean(s$v_ID[ ,x, 1])))
santabarbara_phi <- sapply(1 : dat$N_id, function(x) inv_logit( mean(s$logit_phi) + mean(s$v_ID[ ,x, 2])))
#These values are what we will use to set up our expectations
# For phi, we see a mean of 0.025 and a standard deviation of 0.007
mean(santabarbara_phi)
sd(santabarbara_phi)
range(santabarbara_phi)
# These values are slightly lower than what we had calculated in a previous preregistration (causal) because in there we had used a STAN model where initial attraction for both options was set to zero and birds updated their attractions to both options after each trial. For reference, these were the phi values estimated with that approach
Causal_santabarbara_phi<-c(0.03061692, 0.03430556, 0.04839631, 0.02748937, 0.03125310, 0.03356904, 0.04142868, 0.03397079)
# For lambda, we see a mean of 4.9 and a standard deviation of 2.3
mean(santabarbara_lambda)
sd(santabarbara_lambda)
range(santabarbara_lambda)
# Again, for comparison, here are the lambda values estimated with the prvious approach (they are slightly lower)
Causal_santabarbara_lambda<-c(3.311051, 4.811929, 2.897794, 4.687649, 2.360964, 6.928005, 6.013120, 3.492505)
################################################################################################
# Start the simulation here
########################
# Based on the values for the Santa Barbara birds, we can now set our expectations as inputs for the power simulations. To illustrate the range, a phi of 0.02 means that there is a ~1% shift in attraction if birds would be equally attracted to an option whereas a phi of 0.04 means that there is a ~2% shift.
simulated_sitemean_phi<-c(0.01,0.015,0.02,0.0225,0.025,0.03,0.035,0.04)
simulated_sitesd_phi<-c(0.005,0.005,0.006,0.006,0.007,0.007,0.008,0.008)
simulated_sitemean_lambda<-c(3,4,6,9)
simulated_sitesd_lambda<-c(1,1.5,2.5,3)
#Calculate the number of sites using the number of combinations of the mean and sd options for phi and lambda
totalnumberofsites<-length(simulated_sitemean_phi)*length(simulated_sitemean_lambda)
# We set up a dataframe that will store the relevant information, saving the choice of individuals at each trials and how their attraction scores change correspondingly.
simulatedreversaldata<-matrix(ncol=15,nrow=totalnumberofsites*20*300)
simulatedreversaldata<-as.data.frame(simulatedreversaldata)
colnames(simulatedreversaldata)<-c("Bird_ID", "Session", "Trial", "Reversal","Choice", "CorrectChoice","Phi_mean","Lambda_mean","Site","Phi_sd","Lambda_sd","ThisBirdsPhi","ThisBirdsLambda","Attraction1","Attraction2")
# There are two tubes, lightgrey and darkgrey
# Initially, one of them is the rewarded option (here lightgrey) and the other is the unrewarded option (here darkgrey)
# Rewarded option gives a payoff of 1, unrewarded option gives a payoff of 0
# When a reversal occurs, the rewards are flipped (previous option with 1 now gets 0, previous option with 0 now gets 1)
counter<-1
bird_id<-1
site_id<-1
for (mean_phi in 1:length(simulated_sitemean_phi)) {
current_mean_phi<-simulated_sitemean_phi[mean_phi]
current_sd_phi<-simulated_sitesd_phi[mean_phi]
for (mean_lambda in 1:length(simulated_sitemean_lambda)) {
current_mean_lambda<-simulated_sitemean_lambda[mean_lambda]
current_sd_lambda<-simulated_sitesd_lambda[mean_lambda]
for (bird in 1:20) {
bird_phi<-rnorm(1,current_mean_phi,current_sd_phi)
if(bird_phi<0){bird_phi<-0}
bird_lambda<-rnorm(1,current_mean_lambda,current_sd_lambda)
if(bird_lambda<0){bird_lambda<-0}
# All birds initially have only a small attraction for either option;
# setting these values to zero means that there would be no attraction at all. In this case,
# choosing the wrong option and not getting a reward does not reduce the attraction to that option.
# Accordingly, with values set to zero birds can only learn when they pick the right option.
# Alternatively, setting both these to the same value higher than zero (e.g. 0.1) means
# that individuals initially do not have a preference, but are somewhat attracted to both options.
# In this case, not getting a reward at the wrong option would reduce the attraction to that option.
# Accordingly, with values larger than zero birds can initially potentially learn twice as fast
# because they update attractions both when choosing the correct and the wrong option.
attractionlightgrey<-0.1
attractiondarkgrey<-0.1
# Initially, the light grey option is rewarded
lightgreyreward<-1
darkgreyreward<-0
reversal<-"initial"
# Reflecting the data being collected, the maximum number of sessions a bird can participate in is 30
for (session in 1:30) {
# Each session runs for 10 trials
for (trial in 1:10) {
# We first calculate the probability that the bird will choose the light grey option during this trial based on her relative attraction to this option over the dark grey option
probability_choose_lightgrey<-exp(bird_lambda*attractionlightgrey)/(exp(bird_lambda*attractionlightgrey)+exp(bird_lambda*attractiondarkgrey) )
# Based on the probability, the bird will now make a choice
choselightgrey<-if(runif(1,0,1)<probability_choose_lightgrey){1} else{0}
# For the option they choose, they update their attraction. If they chose the correct option, the learned about the reward and increase their attraction to that option; if they chose the wrong option, they decrease their attraction to that option
if(choselightgrey==1){
attractionlightgrey<-(1-bird_phi)*attractionlightgrey+bird_phi*lightgreyreward
} else {
attractiondarkgrey<-(1-bird_phi)*attractiondarkgrey+bird_phi*darkgreyreward
}
# We store all the relevant information about this trial in the data frame
simulatedreversaldata[counter,]$Bird_ID<-bird_id
simulatedreversaldata[counter,]$Session<-session
simulatedreversaldata[counter,]$Trial<-(session-1)*10+trial
simulatedreversaldata[counter,]$Reversal<-reversal
simulatedreversaldata[counter,]$Site<-site_id
simulatedreversaldata[counter,]$Phi_mean<-current_mean_phi
simulatedreversaldata[counter,]$Lambda_mean<-current_mean_lambda
simulatedreversaldata[counter,]$Phi_sd<-current_sd_phi
simulatedreversaldata[counter,]$Lambda_sd<-current_sd_lambda
simulatedreversaldata[counter,]$Choice<-choselightgrey
simulatedreversaldata[counter,]$Attraction1<-attractionlightgrey
simulatedreversaldata[counter,]$Attraction2<-attractiondarkgrey
simulatedreversaldata[counter,]$ThisBirdsPhi<-bird_phi
simulatedreversaldata[counter,]$ThisBirdsLambda<-bird_lambda
simulatedreversaldata[counter,]$CorrectChoice<-ifelse(reversal=="initial",ifelse(choselightgrey==1,1,0),ifelse(choselightgrey==1,0,1))
counter<-counter+1
#End of trial loop
}
# Birds are done when they have successfully completed the reversal trials after the initial learning trials. Matching the experiments, birds are counted as successful if they chose the correct option in 17 or more of the 20 trials and if they also during those last two sessions chose the correct option 8 or more times out of the 10 times in each session.
if (reversal=="reversal") {
if (sum(simulatedreversaldata[(counter-11):(counter-1),]$CorrectChoice)>7
& sum(simulatedreversaldata[(counter-21):(counter-12),]$CorrectChoice)>7
& sum(simulatedreversaldata[(counter-21):(counter-1),]$CorrectChoice)>16 ) {
break
}
}
# Birds switch from the initial learning to the reversal trials, based on the same definition of being successful. In this case, the reward is now set to be associated with the dark grey option in the following reversal trials.
if(session > 1) {
if ( sum(simulatedreversaldata[(counter-11):(counter-1),]$CorrectChoice)>7
& sum(simulatedreversaldata[(counter-21):(counter-12),]$CorrectChoice)>7
& sum(simulatedreversaldata[(counter-21):(counter-1),]$CorrectChoice)>16 ) {
lightgreyreward<-0
darkgreyreward<-1
reversal<-"reversal"
}
}
#End of session loop
}
bird_id <-bird_id+1
#End of bird loop
}
site_id<-site_id+1
# End of mean lambda loop
}
#End of mean phi loop
}
# Remove unnecessary rows at the bottom that contain only NA values
rem_simulatedreversaldata<-simulatedreversaldata[is.na(simulatedreversaldata$Site)==F,]
# What are the attraction values observed across all individuals in all populations across both initial and reversal for the option that is rewarded during the initial phase
hist(rem_simulatedreversaldata$Attraction1)
# What are the attraction values observed across all individuals in all populations across both initial and reversal for the option that is rewarded during the reversal phase
hist(rem_simulatedreversaldata$Attraction2)
# Boxplots of the 20 individual phi values per population to see whether populations with different mean phi values are different
boxplot(rem_simulatedreversaldata$ThisBirdsPhi~rem_simulatedreversaldata$Phi_mean)
# Save the output - use this code for the modified version with initial attraction scores at 0.1 and individuals only learning about the option they chose
write.csv(rem_simulatedreversaldata,file="SimulatedReversalData_Grackles_PhiLambda_Attraction.csv")
################################################################################################
# End of the simulation
################################################################################################
################################################################################################
# Load previously simulated data
################################################################################################
# These two are the sets we decided on, with initial attractions at 0.1 and eight different phi and four different lambda combinations
simulatedreversaldata_attractionscores_1<-read.csv(url("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/gxpopbehaviorhabitat_SimulatedReversalData_Grackles_PhiLambda_Attraction02_Aug2021.csv"), header=T, sep=",", stringsAsFactors=F)
simulatedreversaldata_attractionscores_2<-read.csv(url("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/gxpopbehaviorhabitat_SimulatedReversalData_Grackles_PhiLambda_Attraction04_Aug2021.csv"), header=T, sep=",", stringsAsFactors=F)
# In both simulations, sites were counted from 1-16; we now assign each site a unique value from 1-32
simulatedreversaldata_attractionscores_2$Site<-simulatedreversaldata_attractionscores_2$Site+16
# We combine the two data sets for the further analyses
simulatedreversaldata_attractionscores<-bind_rows(simulatedreversaldata_attractionscores_1,simulatedreversaldata_attractionscores_2)
################################################################################################
# In the simulations, trials were counted continuously for each bird. We now want to change this so that it restarts counting trials from 1 upward once a bird switches to reversal.
for (birds in 1:length(unique(simulatedreversaldata_attractionscores$Bird_ID))){
currentbird<-unique(simulatedreversaldata_attractionscores$Bird_ID)[birds]
maximuminitial<-max(simulatedreversaldata_attractionscores[simulatedreversaldata_attractionscores$Bird_ID==currentbird & simulatedreversaldata_attractionscores$Reversal == "initial",]$Trial)
simulatedreversaldata_attractionscores[simulatedreversaldata_attractionscores$Bird_ID==currentbird & simulatedreversaldata_attractionscores$Reversal == "reversal",]$Trial<-simulatedreversaldata_attractionscores[simulatedreversaldata_attractionscores$Bird_ID==currentbird & simulatedreversaldata_attractionscores$Reversal == "reversal",]$Trial - maximuminitial
}
# We need to adjust the coding during the reversal learning so that "correct" now matches whether it is correct or not. In the table, choice means which of the two options they chose. If they chose 0, that means they chose the lt gray, if they chose 1, that means they chose the dk gray. CorrectChoice indicates whether this is the rewarded colour or not. It is currently set up so that during the initial learning, option 1 in choice is being rewarded, whereas in the reversal learning, option 0 in choice is being rewarded. To match this better to the original data, we assume that option 1 is rewarded in the intial learning and option 2 is rewarded in the reversal learning (which matches coding as factors 1 and 2).
simulatedreversaldata_attractionscores$Choice <- simulatedreversaldata_attractionscores[simulatedreversaldata_attractionscores$Choice==0,]$Choice<-2
# To use the model to estimate the phi and lambda parameters, we first need to change the column names to match these to the specifications in the model: change Bird_ID to id; change Reversal to Choice, change CorrectChoice to Correct, change Site to Expid
colnames(simulatedreversaldata_attractionscores)<-c("counter","id","Session","Trial","Reversal","Choice","Correct","Phi_mean","Lambda_mean","Site","Phi_sd","Lambda_sd","ThisBirdsPhi","ThisBirdsLambda","Attraction1","Attraction2")
################################################################################################
# Next, we can plot the average choices across trials for the different sites
################################################################################################
# Calculate average performance in reversal learning for different phi values c(0.02,0.0325,0.035,0.45)
se <- function(x) sd(x)/sqrt(length(x))
Prop_correct_reversal_phi_0.02 <- sapply(1:100, function (x) mean(simulatedreversaldata_attractionscores$Correct[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.02]) )
Prop_correct_reversal_SE_phi_0.02 <- sapply(1:100, function (x) se(simulatedreversaldata_attractionscores$Correct[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.02]) )
Prop_correct_reversal_phi_0.02<-1-Prop_correct_reversal_phi_0.02 #reverse the proportion so it is proportion of newly rewarded option correct
Prop_correct_reversal_phi_0.0325 <- sapply(1:100, function (x) mean(simulatedreversaldata_attractionscores$Correct[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.0325]) )
Prop_correct_reversal_SE_phi_0.0325 <- sapply(1:100, function (x) se(simulatedreversaldata_attractionscores$Correct[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.0325]) )
Prop_correct_reversal_phi_0.0325<-1-Prop_correct_reversal_phi_0.0325
Prop_correct_reversal_phi_0.035 <- sapply(1:100, function (x) mean(simulatedreversaldata_attractionscores$Correct[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.035]) )
Prop_correct_reversal_SE_phi_0.035 <- sapply(1:100, function (x) se(simulatedreversaldata_attractionscores$simulatedreversaldata_attractionscores[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.035]) )
Prop_correct_reversal_phi_0.035<-1-Prop_correct_reversal_phi_0.035
Prop_correct_reversal_phi_0.45 <- sapply(1:100, function (x) mean(simulatedreversaldata_attractionscores$Correct[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.45]) )
Prop_correct_reversal_SE_phi_0.45 <- sapply(1:100, function (x) se(simulatedreversaldata_attractionscores$Correct[simulatedreversaldata_attractionscores$Trial==x & simulatedreversaldata_attractionscores$Reversal == "reversal" & simulatedreversaldata_attractionscores$Phi_mean ==0.45]) )
Prop_correct_reversal_phi_0.45<-1-Prop_correct_reversal_phi_0.45
# We now plot the curves that show how the average performance changes with each trial. Initially birds are close to 100% choosing the previously rewarded option. We now want to see how quickly they change from that to never selecting this option. We plot the outcomes with the different values for phi in different colours
plot(Prop_correct_reversal_phi_0.02[1:100],type="n" , lwd=2, lty=1, ylim=c(0,1), ylab="Proportion correct",xlab="Trial number", main = "Simulated Grackle Data")
arrows(1:100,c(Prop_correct_reversal_phi_0.02[1:100])-c(Prop_correct_reversal_SE_phi_0.02[1:100]),1:100,c(Prop_correct_reversal_phi_0.02[1:100], Prop_correct_reversal_phi_0.02[1:100])+c(Prop_correct_reversal_SE_phi_0.02[1:100], Prop_correct_reversal_phi_0.02[1:100]),col="darkgray", code=3, lwd=1, length=0.01, angle = 90)
arrows(1:100,c(Prop_correct_reversal_phi_0.035[1:100])-c(Prop_correct_reversal_SE_phi_0.035[1:100]),1:100,c(Prop_correct_reversal_phi_0.035[1:100], Prop_correct_reversal_phi_0.035[1:100])+c(Prop_correct_reversal_SE_phi_0.035[1:100], Prop_correct_reversal_phi_0.035[1:100]),col="blue", code=3, lwd=1, length=0.01, angle = 90)
#arrows(1:100,c(Prop_correct_reversal_phi_0.45[1:100])-c(Prop_correct_reversal_SE_phi_0.45[1:100]),1:100,c(Prop_correct_reversal_phi_0.45[1:100], Prop_correct_reversal_phi_0.45[1:100])+c(Prop_correct_reversal_SE_phi_0.45[1:100], Prop_correct_reversal_phi_0.45[1:100]),col="red", code=3, lwd=1, length=0.01, angle = 90)
# Now we do the same for lambda
Prop_correct_reversal_lambda_3 <- sapply(1:100, function (x) mean(cleanedsimulatedreversaldata$CorrectChoice[cleanedsimulatedreversaldata$Trial==x & cleanedsimulatedreversaldata$Reversal == "reversal" & cleanedsimulatedreversaldata$Lambda_mean ==3]) )
Prop_correct_reversal_SE_lambda_3 <- sapply(1:100, function (x) se(cleanedsimulatedreversaldata$CorrectChoice[cleanedsimulatedreversaldata$Trial==x & cleanedsimulatedreversaldata$Reversal == "reversal" & cleanedsimulatedreversaldata$Lambda_mean ==3]) )
Prop_correct_reversal_lambda_3<-1-Prop_correct_reversal_lambda_3
Prop_correct_reversal_lambda_45 <- sapply(1:100, function (x) mean(cleanedsimulatedreversaldata$CorrectChoice[cleanedsimulatedreversaldata$Trial==x & cleanedsimulatedreversaldata$Reversal == "reversal" & cleanedsimulatedreversaldata$Lambda_mean ==4.5]) )
Prop_correct_reversal_SE_lambda_45 <- sapply(1:100, function (x) se(cleanedsimulatedreversaldata$CorrectChoice[cleanedsimulatedreversaldata$Trial==x & cleanedsimulatedreversaldata$Reversal == "reversal" & cleanedsimulatedreversaldata$Lambda_mean ==4.5]) )
Prop_correct_reversal_lambda_45<-1-Prop_correct_reversal_lambda_45
# Plot the curves for the different lambda values
plot(Prop_correct_reversal_lambda_3[1:100],type="n" , lwd=2, lty=1, ylim=c(0,1), ylab="Proportion correct",xlab="Trial number", main = "Simulated Grackle Data")
arrows(1:100,c(Prop_correct_reversal_lambda_3[1:100])-c(Prop_correct_reversal_SE_lambda_3[1:100]),1:100,c(Prop_correct_reversal_lambda_3[1:100], Prop_correct_reversal_lambda_3[1:100])+c(Prop_correct_reversal_SE_lambda_3[1:100], Prop_correct_reversal_lambda_3[1:100]),col="darkgray", code=3, lwd=1, length=0.01, angle = 90)
arrows(1:100,c(Prop_correct_reversal_lambda_45[1:100])-c(Prop_correct_reversal_SE_lambda_45[1:100]),1:100,c(Prop_correct_reversal_lambda_45[1:100], Prop_correct_reversal_lambda_45[1:100])+c(Prop_correct_reversal_SE_lambda_45[1:100], Prop_correct_reversal_lambda_45[1:100]),col="red", code=3, lwd=1, length=0.01, angle = 90)
################################################################################################
# Do phi and lambda influence performance?
################################################################################################
performance_mean<-as.data.frame(simulatedreversaldata_attractionscores %>% group_by(id, Phi_mean, Reversal) %>% summarize(mean(Trial)))
boxplot((performance_mean[performance_mean$Reversal=="initial",]$`mean(Trial)`)~performance_mean[performance_mean$Reversal=="initial",]$Phi_mean)
boxplot((performance_mean[performance_mean$Reversal=="reversal",]$`mean(Trial)`)~performance_mean[performance_mean$Reversal=="reversal",]$Phi_mean)
performance_mean_lambda<-as.data.frame(simulatedreversaldata_attractionscores %>% group_by(id, Lambda_mean, Reversal) %>% summarize(mean(Trial)))
boxplot((performance_mean_lambda[performance_mean_lambda$Reversal=="initial",]$`mean(Trial)`)~performance_mean_lambda[performance_mean_lambda$Reversal=="initial",]$Lambda_mean)
boxplot((performance_mean_lambda[performance_mean_lambda$Reversal=="reversal",]$`mean(Trial)`)~performance_mean_lambda[performance_mean_lambda$Reversal=="reversal",]$Lambda_mean)
################################################################################################
# Can we estimate the phi lambda values back from the simulated data
################################################################################################
# We start with a single site, and estimate values for each individual
# Randomly pick one site from the 32 simulated sites
site<-sample(unique(simulatedreversaldata_attractionscores$Site),1)
# Pick the individuals and their trials from that one site
singlesitetest<-simulatedreversaldata_attractionscores[simulatedreversaldata_attractionscores$Site==site ,]
# Renumber the site id so that it is 1 (necessary so that the STAN model knows to work on these data)
singlesitetest$id<-singlesitetest$id-((site-1)*20)
# Change the Choice into whether they picked option 1 or 2
singlesitetest$Choice<-as.integer(as.factor(singlesitetest$Choice))
# If you want to run the estimation only on the reversal trials, use the following, and then in the setting up of the data for the model choose "singlesitetest_reversal" instead of "singlesitetest"
singlesitetest_reversal<-singlesitetest[singlesitetest$Reversal=="reversal",]
# Set up the data for the STAN model
dat_singlesite <- as.list(singlesitetest)
dat_singlesite$N <- nrow(singlesitetest)
dat_singlesite$N_id <- length(unique(singlesitetest$id))
# Define the STAN model. In this model, the initial attraction is set to 0.1 for both options and individuals only update their attraction for the option they chose.
reinforcement_model_nonzeroattraction <- "
data{
int N;
int N_id;
int id[N];
int Trial[N];
int Choice[N];
int Correct[N];
}
parameters{
real logit_phi;
real log_L;
// Varying effects clustered on individual
matrix[2,N_id] z_ID;
vector<lower=0>[2] sigma_ID; //SD of parameters among individuals
cholesky_factor_corr[2] Rho_ID;
}
transformed parameters{
matrix[N_id,2] v_ID; // varying effects on stuff
v_ID = ( diag_pre_multiply( sigma_ID , Rho_ID ) * z_ID )';
}
model{
matrix[N_id,2] A; // attraction matrix
logit_phi ~ normal(0,1);
log_L ~ normal(0,1);
// varying effects
to_vector(z_ID) ~ normal(0,1);
sigma_ID ~ exponential(1);
Rho_ID ~ lkj_corr_cholesky(4);
// initialize attraction scores
for ( i in 1:N_id ) {
A[i,1] = 0.1; A[i,2] = 0.1';
}
// loop over Choices
for ( i in 1:N ) {
vector[2] pay;
vector[2] p;
real L;
real phi;
// first, what is log-prob of observed choice
L = exp(log_L + v_ID[id[i],1]);
p = softmax(L*A[id[i],1:2]' );
Choice[i] ~ categorical( p );
// second, update attractions conditional on observed choice
phi = inv_logit(logit_phi + v_ID[id[i],2]);
pay[1:2] = rep_vector(0,2);
pay[ Choice[i] ] = Correct[i];
A[ id[i] , Choice[i] ] = ( (1-phi)*(A[ id[i] , Choice[i] ]) + phi*pay[Choice[i]])';
}//i
}
"
# Now we run the STAN model on the data from this one site
m_singlesite <- stan( model_code = reinforcement_model_nonzeroattraction, data=dat_singlesite ,iter = 5000, cores = 4, chains=4, control = list(adapt_delta=0.9, max_treedepth = 12))
# We extract samples from the posterior of the STANmodel
posterior_singlesite<-extract.samples(m_singlesite)
# We calculate what lambda and phi values the STAN model estimated for each individual from this one site
lambda_singlesite <- sapply(1 : dat_singlesite$N_id, function(x) exp( mean(posterior_singlesite$log_L) + mean(posterior_singlesite$v_ID[ ,x, 1])))
phi_singlesite <- sapply(1 : dat_singlesite$N_id, function(x) inv_logit( mean(posterior_singlesite$logit_phi) + mean(posterior_singlesite$v_ID[ ,x, 2])))
# We set up a quick projection function to see how well the estimated STAN values can recover the performance of the birds at this one simulated site
Projection_fct <- function(Tmax = 100, N_id = 20){
#Create output object
d_Overall <- c()
#Loop over birds
for (ind in 1:N_id) {
#Set initial attractions to 0
A <- c(0.1,0.1)
#Create output matrix to record choices and payoffs
d <- data.frame(id=ind, trial = 1:Tmax, Choice=NA, Payoff=NA, Experiment = NA)
d$Experiment[1:50] <- 0
d$Experiment[51:100] <- 1
# Start simulation loop
for (i in 1:Tmax) {
Prob <- c()
for (xx in 1:2) {
Prob[xx] <- exp(lambda_singlesite[ind]*A[xx]) / sum(exp(lambda_singlesite[ind]*A))
}
#Make choice proportional to attraction scores
d$Choice[which(d$id==ind & d$trial==i)] <- sample(c(1:2), size = 1, prob = Prob)
pay <- c(0,0)
if (d$Experiment[i] == 0){
if (d$Choice[i] == 1){
d$Payoff[i] <- 1
pay[1] <- 1
} else {
d$Payoff[i] <- 0
}
} else {
if (d$Choice[i] == 1){
d$Payoff[i] <- 0
} else {
d$Payoff[i] <- 1
pay[2] <- 1
}
}
#Update Attractions
for (xx in 1:2) {
A[xx] <- (1-phi_singlesite[ind]) * A[xx] + phi_singlesite[ind] * pay[xx]
}
}#time i
d_Overall <- rbind(d_Overall, d)
}#ind
return(d_Overall)
}
# We run this projection to get the expected performance of the birds based on the STAN estimates
dat_sim <- Projection_fct()
# We now plot the learning curves for this birds at this one simulated site and the projected learning curves based on the STAN estimates
#Compute mean correct and SE of choices for simulated and projected data
Prop_correct_sim<- c()
se <- function(x) sd(x)/sqrt(length(x))
Prop_correct_sim <- sapply(1:100, function (x) mean(dat$Payoff[dat$trial==x]) )
Prop_correct_sim_SE <- sapply(1:100, function (x) se(dat$Payoff[dat$trial==x]) )
Prop_correct_initial<- c()
Prop_correct_initial <- sapply(1:100, function (x) mean(singlesitetest$Correct[singlesitetest$Trial==x & singlesitetest$Reversal == "initial"]) )
Prop_correct_initial_SE <- sapply(1:100, function (x) se(singlesitetest$Correct[singlesitetest$Trial==x & singlesitetest$Reversal == "initial"]) )
Prop_correct_reversal<- c()
Prop_correct_reversal <- sapply(1:100, function (x) mean(singlesitetest$Correct[singlesitetest$Trial==x & singlesitetest$Reversal == "reversal"]) )
Prop_correct_reversal_SE <- sapply(1:100, function (x) se(singlesitetest$Correct[singlesitetest$Trial==x & singlesitetest$Reversal == "reversal"]) )
op<-par()
par(mfrow=c(2,1),
oma=c(1.1,2.5,0,0),
mar=c(2.5,2.4,1.5,0.1))
plot(c(Prop_correct_initial[1:50], Prop_correct_reversal[1:50]), type="n", lwd=2, lty=1, ylim=c(0,1), ylab="",xlab="", main = "simulated site 40")
arrows(1:100,c(Prop_correct_initial[1:50], Prop_correct_reversal[1:50])-c(Prop_correct_initial_SE[1:50], Prop_correct_reversal_SE[1:50]),1:100,c(Prop_correct_initial[1:50], Prop_correct_reversal[1:50])+c(Prop_correct_initial_SE[1:50], Prop_correct_reversal_SE[1:50]),col="darkgray", code=3, lwd=1, length=0.01, angle = 90)
lines(c(Prop_correct_initial[1:50], Prop_correct_reversal[1:50]), lwd=2)
abline(v = 51, lty =2)
plot(Prop_correct_sim, type="n", lwd=2, lty=1, ylim=c(0,1), ylab="",xlab="", main = "Simulated Data")
arrows(1:100,Prop_correct_sim-Prop_correct_sim_SE,1:100,Prop_correct_sim+Prop_correct_sim_SE, code=3, lwd=1,col="darkgray", length=0.01, angle = 90)
lines(Prop_correct_sim, lwd=2)
abline(v = 51, lty =2)
mtext(side = 1, "Trial", line = 0, outer = TRUE, cex = 1.2)
mtext(side = 2, "Proportion Correct", line = 0.5, outer = TRUE, cex = 1.2)
par(op)
########################################################################################
# Now we run a STAN model with all simulated sites, estimating both means of phi and lambda for each site as well as each individuals' phi and lambda
# We set up the data for the model
dat_full <- as.list(simulatedreversaldata_attractionscores)
dat_full$N <- nrow(simulatedreversaldata_attractionscores)
dat_full$N_id <- length(unique(simulatedreversaldata_attractionscores$id))
dat_full$N_exp <- length(unique(simulatedreversaldata_attractionscores$Site))
dat_full$Choice <- as.numeric(as.factor(dat_full$Choice))
# This STAN model, in addition to estimating phi and lambda for each individual, also estimates means for each site. It again starts with attractions set to 0.1 and assumes that individuals only learn about the option they chose.
reinforcement_model_id_site_nonzeroattraction <- "
data{
int N;
int N_id;
int N_exp;
int id[N];
int Site[N];
int Trial[N];
int Choice[N];
int Correct[N];
}
parameters{
real logit_phi;
real log_L;
// Varying effects clustered on individual
matrix[2,N_id] z_ID;
vector<lower=0>[2] sigma_ID; //SD of parameters among individuals
cholesky_factor_corr[2] Rho_ID;
// Varying effects clustered on experimenter
matrix[2,N_exp] z_EXP;
vector<lower=0>[2] sigma_EXP;
cholesky_factor_corr[2] Rho_EXP;
}
transformed parameters{
matrix[N_id,2] v_ID; // varying effects on individuals
matrix[N_exp,2] v_EXP; // varying effects on experimenter
v_ID = ( diag_pre_multiply( sigma_ID , Rho_ID ) * z_ID )';
v_EXP = ( diag_pre_multiply( sigma_EXP , Rho_EXP ) * z_EXP )';
}
model{
matrix[N_id,2] A; // attraction matrix
logit_phi ~ normal(0,1);
log_L ~ normal(0,1);
// varying effects
to_vector(z_ID) ~ normal(0,1);
sigma_ID ~ exponential(1);
Rho_ID ~ lkj_corr_cholesky(4);
to_vector(z_EXP) ~ normal(0,1);
sigma_EXP ~ exponential(1);
Rho_EXP ~ lkj_corr_cholesky(4);
// initialize attraction scores
for ( i in 1:N_id ) {
A[i,1] = 0.1; A[i,2] = 0.1';
}
// loop over Choices
for ( i in 1:N ) {
vector[2] pay;
vector[2] p;
real L;
real phi;
// first, what is log-prob of observed choice
L = exp(log_L + v_ID[id[i],1]+ v_EXP[Site[i],1]);
p = softmax(L*A[id[i],1:2]' );
Choice[i] ~ categorical( p );
// second, update attractions conditional on observed choice
phi = inv_logit(logit_phi + v_ID[id[i],2]+ v_EXP[Site[i],2]);
pay[1:2] = rep_vector(0,2);
pay[ Choice[i] ] = Correct[i];
A[ id[i] , Choice[i] ] = ( (1-phi)*(A[ id[i] , Choice[i] ]) + phi*pay[Choice[i]])';
}//i
}
"
# Now we run the STAN model on the full data. This will take a few hours.
m_simulated_full <- stan( model_code = reinforcement_model_id_site_nonzeroattraction, data=dat_full ,iter = 5000, cores = 4, chains=4, control = list(adapt_delta=0.9, max_treedepth = 12))
# We extract the posterior from the STAN model
s <- extract.samples(m_simulated_full)
# We calculate the mean phi and lambda for each site as estimated by the STAN model
lambda <- sapply(1 : dat_full$N_exp, function(x) exp( mean(s$log_L) + mean(s$v_EXP[ ,x, 1])))
phi <- sapply(1 : dat_full$N_exp, function(x) inv_logit( mean(s$logit_phi) + mean(s$v_EXP[ ,x, 2])))
# We calculate the original mean phi and lambda values for each site. Because we randomly generated 20 individuals per site, these original values will be slightly different than the means we set as starting values for the simulations
originallambda<-as.vector(as.data.frame(simulatedreversaldata_attractionscores %>% group_by(Site) %>% summarise(mean(ThisBirdsLambda)))[,2])
originalphi<-(as.data.frame(simulatedreversaldata_attractionscores %>% group_by(Site) %>% summarise(mean(ThisBirdsPhi)))[,2])
# We want to know whether the STAN model estimated two sites to be different or not. As criterion for this, we check whether the samples in the posterior for pairs of site overlap or not. That means that, when calculating the difference between the samples in the posterior for a pair of site, they are assumed to be different if the differences are all on one side of zero.
# We first check our power to estimate difference in mean phi between sites
# We calculate the pairwise differences between sites for the estimated phi means
pairwisedifferences_phiposteriors<-list(sapply(1 : (dat_full$N_exp * (dat_full$N_exp-1)/2), function(x) (s$v_EXP[ ,combn(1:dat_full$N_exp,m=2)[1,x], 2]-(s$v_EXP[ ,combn(1:dat_full$N_exp,m=2)[2,x], 2]))))
# We convert it into a dataframe for easier manipulation
estimates_pairwisedifferences_phiposteriors<-as.data.frame(precis(pairwisedifferences_phiposteriors,depth=2))
# We add a column that classifies as true or false whether the estimates crosses zero (as opposed to all values being either smaller or larger than zero)
estimates_pairwisedifferences_phiposteriors$crosseszero<-ifelse(estimates_pairwisedifferences_phiposteriors$mean>0,estimates_pairwisedifferences_phiposteriors$`5.5%`<0,estimates_pairwisedifferences_phiposteriors$`94.5%`>0)
# We go back and again calculate the original phi values for each site (the average across the random sample of 20 birds simulated at each site)
eachsimulatedsitesobservedphi<-as.data.frame(simulatedreversaldata_attractionscores %>% group_by(Site) %>% summarise(mean(ThisBirdsPhi)))
colnames(eachsimulatedsitesobservedphi)<-c("Site","Site_ObservedPhi_mean")
# We want to know all pairwise differences in the original means across the 32 sites
observedphi_combinations_wide<-as.data.frame((combn(eachsimulatedsitesobservedphi$Site_ObservedPhi_mean,m=2)))
# observedphi_combinations_wide is a long vector with two rows. To add it to the data frame we created before that has the difference in the estimated means from the STAN model, we need to transverse it into a format with two columns
observedphi_combinations_long<-matrix(ncol=2,nrow=(dat_full$N_exp * (dat_full$N_exp-1)/2))
observedphi_combinations_long<-as.data.frame(observedphi_combinations_long)
observedphi_combinations_long[,1]<-as.numeric(observedphi_combinations_wide[1,])
observedphi_combinations_long[,2]<-as.numeric(observedphi_combinations_wide[2,])
colnames(observedphi_combinations_long)<-c("observedphi_firstsiteforcomparison","observedphi_secondsiteforcomparison")
observedphi_combinations_long$observedsitedifference<-abs(observedphi_combinations_long$observedphi_firstsiteforcomparison-observedphi_combinations_long$observedphi_secondsiteforcomparison)
# We can now add the column from observedphi_combinations_long that has the differences in the phi means in the original data to the dataframe that has the differences from the STAN estimations
estimates_pairwisedifferences_phiposteriors$observedsitedifference<-observedphi_combinations_long$observedsitedifference
# As a first check, we can see whether the cases where the STAN model estimates that the two sites are not different (difference crosses zero is TRUE) show smaller differences in the original phis of the two sites than those cases where the STAN model does estimate the difference to be clear.
boxplot(estimates_pairwisedifferences_phiposteriors$simulatedsitedifference~estimates_pairwisedifferences_phiposteriors$crosseszero)
# We now also want to add the information about how different two sites were supposed to be based on the input values of the phi means in the simulations. We assumed variation at each site, so when we randomly simulated 20 individuals the mean across their phis can be lower or higher than what we set as the input. We can see that in the following plot, where on the x-axis we have the differences between sites that were set in the inputs and on the y-axis the actual differences in the original phi means that we obtained in the simulations.
plot(dat_phidifferences$differencedetected~combined_estimates_pairwisedifferences_phiposteriors$observedsitedifference)
# We first obtain the input value for each site
eachsimulatedsitessimulatedphi<-as.data.frame(simulatedreversaldata_attractionscores %>% group_by(Site) %>% summarise(mean(Phi_mean)))
colnames(eachsimulatedsitessimulatedphi)<-c("Site","Site_SimulatedPhi_mean")
# We set up all pairwise comparisons between sites again
simulatedphi_combinations_wide<-as.data.frame((combn(eachsimulatedsitessimulatedphi$Site_SimulatedPhi_mean,m=2)))
# The simulatedphi_combinations_wide needs to again be transformed into a data frame with the comparisons in columns
simulatedphi_combinations_long<-matrix(ncol=2,nrow=(dat_full$N_exp * (dat_full$N_exp-1)/2))
simulatedphi_combinations_long<-as.data.frame(simulatedphi_combinations_long)
simulatedphi_combinations_long[,1]<-as.numeric(simulatedphi_combinations_wide[1,])
simulatedphi_combinations_long[,2]<-as.numeric(simulatedphi_combinations_wide[2,])
colnames(simulatedphi_combinations_long)<-c("simulatedphi_firstsiteforcomparison","simulatedphi_secondsiteforcomparison")
simulatedphi_combinations_long$simulatedsitedifference<-abs(simulatedphi_combinations_long$simulatedphi_firstsiteforcomparison-simulatedphi_combinations_long$simulatedphi_secondsiteforcomparison)
# We add the column from simulatedphi_combinations_long that has the differences in the preset phi means for each simulated site to the dataframe that has the differences from the STAN estimations
estimates_pairwisedifferences_phiposteriors$simulatedsitedifference<-simulatedphi_combinations_long$simulatedsitedifference
# We now have a dataframe with all the relevant information to determine our power to detect differences between sites when we have a sample of 20 individuals from each site. To estimate power, we run additional STAN models that estimate how the probability that a difference between a pair of sites crosses zero (meaning that we cannot tell the sites apart) depends on the difference in the phi means between those sites.
# We set up the data for the STAN model. There are two predictors: first, the differences in the inputs for the simulations (converted to categories) and second the differences in the original phi means at each of the simulated sites (log converted because we assume that the probability to detect a difference will drop off exponentially rather than change linearly)
dat_phidifferences <- list(differencedetected = as.integer(estimates_pairwisedifferences_phiposteriors$crosseszero),
simulatedsitedifference = as.integer(as.factor( estimates_pairwisedifferences_phiposteriors$simulatedsitedifference)),
observedsitedifference = standardize(estimates_pairwisedifferences_phiposteriors$observedsitedifference)
)
dat_phidifferences$inv_observedsitedifference<-log(estimates_pairwisedifferences_phiposteriors$observedsitedifference)
# We run the first model that estimates how the probability to detect a difference depends on differences in the input values for the simulated sites
m_phidifferences <- ulam(
alist(
differencedetected~dbinom(1,p),
logit(p) <- b[simulatedsitedifference],
b[simulatedsitedifference]~dnorm(0,1.5)
),data=dat_phidifferences, log_lik=TRUE, messages=FALSE,cmdstan=T)
precis( m_phidifferences , depth=2 )
# We can plot the estimated probabilities to detect a difference across the scale of differences in the input values. The input values were set up such that they can be apart by one of 9 differences
plot(inv_logit(precis(m_phidifferences,depth=2)[,1])~c(0,0.0025,0.005,0.0075,0.01,0.015,0.02,0.025,0.03), lwd=2, lty=1, ylim=c(0,1),xlim=c(0,0.04), ylab=strwrap("Probability model estimates difference could be zero"),xlab="Difference of mean phis between pairs of sites", main = "Power to detect site differences with n=20 birds/site",frame.plot=F)
arrows(x0=c(0,0.0025,0.005,0.0075,0.01,0.015,0.02,0.025,0.03),x1=c(0,0.0025,0.005,0.0075,0.01,0.015,0.02,0.025,0.03),y0=inv_logit(precis(m_phidifferences,depth=2)[,3]),y1=inv_logit(precis(m_phidifferences,depth=2)[,4])
,col=col.alpha(rangi2,0.3), code=3, lwd=5, length=0.01, angle = 90)
# Next, we run the second model that estimates how the probability to detect a difference depends on differences in the original mean values across the 20 random birds at each of the simulated sites
m_phidifferences_inv <- ulam(
alist(
differencedetected~dbinom(1,p),
logit(p) <- a+b*inv_observedsitedifference,
b~dnorm(0,1.5),
a~dnorm(0,1.5)
),data=dat_phidifferences, log_lik=TRUE, messages=FALSE,cmdstan=T)
# To plot the change in probability across the continuous differences in the original mean values for the sites, we set up a dataframe that has the values with the differences in the original mean values and use the model to calculate the probability for each of these values
sequence_mean_differences<-c(sort(unique(dat_phidifferences$inv_observedsitedifference)))
sequence_mean_differences
slopeestimate<-c()
for (i in 1:length(sequence_mean_differences)) {
slopeestimate[i]<-inv_logit(precis( m_phidifferences_inv , depth=2 )[2,1]+precis( m_phidifferences_inv , depth=2 )[1,1]*sequence_mean_differences[i])
}
inv_sequence_mean_differences<-exp(sequence_mean_differences)
plot(slopeestimate~inv_sequence_mean_differences, lwd=2, lty=1, ylim=c(0,1),xlim=c(0,0.04), ylab=strwrap("Probability model estimates difference could be zero"),xlab="Difference of mean phis between pairs of sites", main = "Power to detect site differences with n=20 birds/site",frame.plot=F)
mu <- link( m_phidifferences_inv , data=data.frame( inv_observedsitedifference=sequence_mean_differences ) )
mu_mean <- apply( mu , 2 , mean )
mu_ci <- apply( mu , 2 , PI , prob=0.97 )
lines( inv_sequence_mean_differences , mu_mean , lwd=2 )
shade( mu_ci , inv_sequence_mean_differences , col=col.alpha(rangi2,0.3) )
#####################################################################################################
# We next check our power to estimate difference in mean lambda between sites
# We calculate the pairwise differences between sites for the estimated lambda means
pairwisedifferences_lambdaposteriors<-list(sapply(1 : (dat_full$N_exp * (dat_full$N_exp-1)/2), function(x) (s$v_EXP[ ,combn(1:dat_full$N_exp,m=2)[1,x], 1]-(s$v_EXP[ ,combn(1:dat_full$N_exp,m=2)[2,x], 1]))))
# We convert it into a dataframe for easier manipulation
estimates_pairwisedifferences_lambdaposteriors<-as.data.frame(precis(pairwisedifferences_lambdaposteriors,depth=2))
# We add a column that classifies as true or false whether the estimates crosses zero (as opposed to all values being either smaller or larger than zero)
estimates_pairwisedifferences_lambdaposteriors$crosseszero<-ifelse(estimates_pairwisedifferences_lambdaposteriors$mean>0,estimates_pairwisedifferences_lambdaposteriors$`5.5%`<0,estimates_pairwisedifferences_lambdaposteriors$`94.5%`>0)
# We go back and again calculate the original phi values for each site (the average across the random sample of 20 birds simulated at each site)
eachsimulatedsitesobservedlambda<-as.data.frame(simulatedreversaldata_attractionscores %>% group_by(Site) %>% summarise(mean(ThisBirdsLambda)))