Skip to content
Permalink
master
Switch branches/tags
Go to file
 
 
Cannot retrieve contributors at this time
---
title: "Is behavioral flexibility linked with exploration, but not boldness, persistence, or motor diversity?"
author:
- 'McCune KB^1^'
- 'MacPherson M^1^'
- 'Rowney C^2^'
- 'Bergeron L^1^'
- 'Folsom M^2^'
- 'Deffner D^2^'
- 'Breen A^2^'
- '[Logan CJ](http://CorinaLogan.com)^2^'
date: '`r Sys.Date()`'
output:
html_document:
toc: true
toc_depth: 4
toc_float:
collapsed: false
code_folding: hide
md_document:
toc: true
pdf_document:
keep_tex: yes
latex_engine: xelatex
github_document:
toc: true
bibliography: MyLibrary.bib
header-includes:
- \usepackage[left]{lineno}
- \linenumbers
---
##### Affiliations:
1) University of California Santa Barbara
2) Max Planck Institute for Evolutionary Anthropology
*Corresponding author: kelseybmccune@gmail.com
```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=TRUE)
knitr::opts_chunk$set(echo = TRUE)
```
<img width="50%" src="logoPCIecology.png">
**Cite as:** McCune K, MacPherson M, Rowney C, Bergeron L, Folsom M, Logan CJ. 2019. [Is behavioral flexibility linked with exploration, but not boldness, persistence, or motor diversity?](http://corinalogan.com/Preregistrations/g_exploration.html) (http://corinalogan.com/Preregistrations/g_exploration.html) In principle acceptance by *PCI Ecology* of the version on 27 Mar 2019 [https://github.com/corinalogan/grackles/blob/master/Files/Preregistrations/g_exploration.Rmd](https://github.com/corinalogan/grackles/blob/master/Files/Preregistrations/g_exploration.Rmd).
<img width="5%" src="logoOpenAccess.png"> <img width="5%" src="logoOpenCode.png"> <img width="5%" src="logoOpenPeerReview.png">
**This preregistration has been pre-study peer reviewed and received an In Principle Recommendation by:**
Jeremy Van Cleve (2019) Probing behaviors correlated with behavioral flexibility. *Peer Community in Ecology*, 100020. [10.24072/pci.ecology.100020](https://ecology.peercommunityin.org/public/rec?id=29&reviews=True)
- Reviewers: two anonymous reviewers
\newpage
## ABSTRACT
This is one of the first studies planned for our long-term research on the role of behavioral flexibility in rapid geographic range expansions. **Project background:** Behavioral flexibility, the ability to change behavior when circumstances change based on learning from previous experience (@mikhalevich_is_2017), is thought to play an important role in a species' ability to successfully adapt to new environments and expand its geographic range (e.g., [@lefebvre1997feeding], [@griffin2014innovation], [@chow2016practice], [@sol2000behavioural], [@sol2002behavioural], [@sol2005big]). However, behavioral flexibility is rarely directly tested at the individual level, thus limiting our ability to determine how it relates to other traits, which limits the power of predictions about a species' ability to adapt behavior to new environments. We use great-tailed grackles (a bird species) as a model to investigate this question because they have rapidly expanded their range into North America over the past 140 years ([@wehtje2003range], [@peer2011invasion]) (see an overview of the [5-year project timeline](https://github.com/corinalogan/grackles/blob/master/README.md)). **This investigation**: In this piece of the long-term project, we aim to understand whether grackle behavioral flexibility (color tube reversal learning - described in a separate [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md)) correlates (or not) with individual differences in the exploration of new environments and novel objects, boldness, persistence, and motor diversity (and whether the flexibility manipulation made such correlations more detectable). Results will indicate whether consistent individual differences in these traits might interact with measures of flexibility (reversal learning and solution switching). This will improve our understanding of which variables are linked with flexibility and how they are related, thus putting us in an excellent position to further investigate the mechanisms behind these links in future research.
### [Video summary](https://youtu.be/Xd_nYV9Lj7E) [https://youtu.be/Xd_nYV9Lj7E](https://youtu.be/Xd_nYV9Lj7E)
## INTRODUCTION
## ASSOCIATED PREREGISTRATION
Our hypotheses, methods, and analysis plans are described in the peer-reviewed preregistration of this article, which is included below as the [Methods](#methods).
## DEVIATIONS FROM THE PREREGISTRATION
**After data collection began and before data analysis:**
1) We added interobserver reliability analyses for the response variables to determine how repeatable our data collection methods were.
**During data analysis:**
2) Models using latency as a response variable were changed from Poisson to Gaussian because latency is a continuous variable
3) For the repeatability analyses, we additionally analyzed the data using the rptR package to ensure results were consistent with the MCMCglmm results.
## RESULTS
## DISCUSSION
## METHODS
Below is the preregistration that passed pre-study peer review.
### A. STATE OF THE DATA
NOTE: all parts of the preregistration are included in this one manuscript.
**Prior to collecting any data:** This preregistration was written and submitted to PCI Ecology for peer review (Sep 2018).
**After data collection had begun (and before any data analysis was conducted):** This preregistration was peer reviewed at PCI Ecology, revised, and resubmitted (Feb 2019), and passed pre-study peer review (Mar 2019). See the [peer review history](https://ecology.peercommunityin.org/public/rec?id=29&reviews=True). Interobserver reliability analyses were added (Feb 2021).
### B. PARTITIONING THE RESULTS
We may decide to present the results from different hypotheses in separate papers.
### C. HYPOTHESES
#### H1: [Behavioral flexibility](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md) (indicated by individuals that are faster at functionally changing their behavior when circumstances change; measured by reversal learning and switching between options on a multi-access box) is positively correlated with the exploration of new environments and novel objects, but not with other behaviors (i.e., boldness, persistence, or motor diversity) (see @mikhalevich_is_2017 for theoretical background about our flexibility definition).
We will first verify that our measures of exploration, boldness and persistence represent repeatable, inherent individual differences in behavior (i.e., personality). Individuals show consistent individual differences in behavior if the variance in latency to approach the task is smaller within individuals compared to variance in latency among individuals (for exploration and boldness assays). The same definition applies to persistence with the number of touches as the measured variable. If there is no repeatibility of these behaviors within individuals, then performance is likely state dependent (e.g., it depends on their fluctuating motivation, hunger levels, etc.) and/or reliant on the current context of the tasks.
**Predictions 1-5:** Individuals in the experimental group where flexibility (as measured by reversal learning and on a multi-access box) was manipulated (such that individuals in the manipulated group became faster at switching) will be more exploratory of new environments (P1; methods similar to free-entry open field test as in @mettke2009spatial) and novel objects (P2; methods as in @mettke2009spatial) than individuals in the control group where flexibility was not increased, and there will be no difference between the groups in persistence (P3), boldness (P4; methods as in @logan2016behavioral), or motor diversity (P5) (as found in @logan2016behavioral). We do not expect the flexibility manipulation to causally change the nature of the relationship between flexibility and any of the other measured variables. Instead, we expect the manipulation to potentially enhance individual variation, thus making it easier for us to detect a correlation if one exists.
**P1-P5 alternative:** If the flexibility manipulation does not work in that those individuals in the experimental condition are not more flexible than control individuals, then we will analyze the individuals from both conditions as one group. In this case, we will assume that we were not able to influence their flexibility and that whatever level of flexibility they had coming into the experiment reflects the general individual variation in the population. This experiment will then elucidate whether general individual variation in flexibility relates to exploratory behaviors. The predictions are the same as above. The following alternatives apply to both cases: if the manipulation works (in which case we expect stronger effects for the manipulated group), and if the manipulation doesn't work (in which case we expect individuals to vary across all of the measured variables and for these variables to potentially interact).
**P1 alternative 1:** There is a positive correlation between exploration and both dependent variables in reversal learning (one accounts for exploration in reversal learning [the ratio] and the other does not). This suggests that flexibility is not independent of exploration and could indicate that another trait is present that could be explaining individual variation in flexibility as well as in exploration. This other trait or traits could be something such as boldness or persistence.
**P1 alternative 2a:** There is a positive correlation between exploration and the dependent variable that does not account for exploration (number of trials to reverse), but not the flexibility ratio, which suggests that performance overall in reversal learning is partially explained by variation in exploration, but that flexibility and exploration are separate traits because using a measure that accounts for exploration still shows variation in flexibility.
**P1 alternative 2b:** There is a negative correlation between exploration and the flexibility ratio that accounts for exploration, but not with the number of trials to reverse. This could be an artifact of accounting for exploration in both variables.
**P1 alternative 3:** There is no correlation between exploration and either dependent variable in reversal learning. This indicates that both dependent variables measure traits that are independent of exploration.
**P1 alternative 4:** There is no correlation between exploration and either dependent variable in reversal learning because our novel object and novel environment methods are inappropriate for measuring exploratory tendency. These measures of exploration both incorporate novelty and thus may measure boldness rather than exploration. This is supported by a positive correlation between behavioral responses to our exploration and boldness assays.
**P3 alternative 1:** There is a positive correlation between persistence and the number of incorrect choices in reversal learning before making the first correct choice. This indicates that individuals that are persistent in one context are also persistent in another context.
**P3 alternative 2:** There is no correlation between persistence and the number of incorrect choices in reversal learning before making the first correct choice. This indicates that flexibility is an independent trait.
![Figure 1.](g_exploreFig1.png)
***Figure 1.*** An overview of the study design and a selection of the variables we will measure for each assay. Exploration will be measured by comparing individual behavior within a familiar environment to behavior towards a novel environment, as well as response to a familiar object vs. a novel object within the familiar environment that contains their regular food. Boldness will be measured as the willingness to eat next to a threatening object (familiar, novel oject, or a taxidermic predator) in their familiar environment. Persistence will be measured as the number of touches to the novel environment and novel object in the Exploration assay, the objects in the Boldness assay, and the multi-access box in a separate [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md). Motor diversity will be measured using the multi-access box in a separate [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md). After the flexibility manipulation occurs, assays will be conducted at least twice (e.g., Time 1, Time 2) and differences (if any) between the control and manipulated groups in the behavioral flexibility [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md) will be compared across time and, with persistence, across tests (e.g., Test 1, Test 2) because persistence is measured in four different assays.
#### H2: Captive and wild individuals may respond differently to assays measuring exploration and boldness.
**P6:** Individuals assayed while in captivity are less exploratory and bold than when they are again assayed in the wild, and as compared to separate individuals assayed in the wild, potentially because captivity is an unfamiliar situation.
**P6 alternative 1:** Individuals in captivity are more exploratory and bold than wild individuals (testing sessions matched for season), and captive individuals show more exploratory and bold behaviors than when they are subsequently tested in the wild, potentially because the captive environment decreases the influence of predation, social interactions and competition.
**P6 alternative 2:** There is no difference in exploration and boldness between individuals in captivity and individuals in the wild (matched for season), potentially because in both contexts our data is biased by sampling only the types of individuals that were most likely to get caught in traps.
**P6 alternative 3:** Captive individuals, when tested again after being released, show no difference in exploratory and bold behaviors because our methods assess inherent personality traits that are consistent across the captive and wild contexts in this taxa.
### D. METHODS
**Planned Sample**
Great-tailed grackles are caught in the wild in Tempe, Arizona USA for individual identification (colored leg bands in unique combinations). Some individuals (~32) are brought temporarily into aviaries for testing, and then they will be released back to the wild. Grackles are individually housed in an aviary (each 244cm long by 122cm wide by 213cm tall) at Arizona State University for a maximum of three months where they have ad lib access to water at all times and are fed Mazuri Small Bird maintenance diet ad lib during non-testing hours (minimum 20h per day), and various other food items (e.g., peanuts, grapes, bread) during testing (up to 3h per day per bird). Individuals are given three to four days to habituate to the aviaries and then their test battery begins on the fourth or fifth day (birds are usually tested six days per week, therefore if their fourth day in the aviaries occurs on a day off, then they are tested on the fifth day instead). For hypothesis 2 we will attempt to test all grackles in the wild that are color-banded.
**Sample size rationale**
We will test as many birds as we can in the approximately three years at this field site given that the birds only participate in tests in aviaries during the non-breeding season (approximately September through March). The minimum sample size for captive subjects will be 16, however we expect to be able to test up to 32 grackles in captivity. We catch grackles with a variety of methods, some of which decrease the likelihood of a selection bias for exploratory and bold individuals because grackles cannot see the traps (i.e. mist nets). In samplng all banded birds in the wild, we will therefore have a better idea of the variation in exploration and boldness behaviors in this population.
**Data collection stopping rule**
We will stop testing birds once we have completed two full aviary seasons (likely in March 2020) if the sample size is above the minimum suggested boundary based on model simulations (see section "[Ability to detect actual effects](#Ability-to-detect-actual-effects)" below). If the minimum sample size is not met by this point, we will continue testing birds at our next field site (which we move to in the summer of 2020) until we meet the minimum sample size.
#### **Open materials**
[Testing protocols](https://docs.google.com/document/d/1sEMc5z2fw6S9C-wVfc2zV331CRPpu3NuA7IhSFUZJpE/edit?usp=sharing) for exploration of new environments and objects, boldness, persistence, and motor diversity.
#### **Open data**
When the study is complete, the data will be published in the Knowledge Network for Biocomplexity's data repository.
#### **Randomization and counterbalancing**
There is no randomizing. The order of the three tasks will be counterbalanced across birds (using https://www.random.org to randomly assign individuals to one of three experimental orders).
1/3 of the individuals will experience:
1. Exploration environment
2. Exploration object
3. Boldness
1/3 of the individuals will experience:
1. Exploration object
2. Boldness
3. Exploration environment
1/3 of the individuals will experience:
1. Boldness
2. Exploration environment
3. Exploration object
#### **Blinding of conditions during analysis**
No blinding is involved in this study. NOTE Feb 2021: interobserver reliability analyses were conducted with hypothesis-blind video coders.
#### **Variables included in analyses 1-5**
NOTE: to view a list of these variables in a table format, please see our Google [sheet](https://docs.google.com/spreadsheets/d/1nhFkqTFWeAeWli8FU8n7mDiWGBuCeduzf8tWN3wPQeE/edit?usp=sharing), which describes whether they are a dependent variable (DV), independent variable (IV), or random effect (RE). Note: when there is more than one DV per model, all models will be run once per DV.
**ANALYSIS 1** - *REPEATABILITY of boldness, persistence and exploration*
**Dependent variables**
1) Boldness: Latency to land on the table - OR - Latency to eat the food - OR - Latency to touch a threatening object next to food (we will choose the variable with the most data)
2) Persistence: Number of touches to an apparatus per time (multi-access box in the behavioral flexibility [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md), novel environment in P1, and objects in P2 and P4)
3) Exploration of novel environment: Latency to enter a novel environment set inside a familiar environment
4) Exploration of novel object: Latency to land on the table next to an object (novel, familiar) (that does not contain food) in a familiar environment (that contains maintenance diet away from the object) - OR - latency to touch an object (novel, familiar) (choose the variable with the most data)
**Independent variables**
1) [Condition](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md): control, flexibility manipulation
2) ID (random effect because multiple measures per individual)
**ANALYSIS 2** - *H1: P1-P5: flexibility correlates with exploratory behaviors*
**Dependent variables**
1) The **number of trials to reverse** a preference in the last reversal that individual participated in (an individual is considered to have a preference if it chose the rewarded option at least 17 out of 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)). See behavioral flexibility [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md) for details.
2) If the number of trials to reverse a preference does not positively correlate with the number of trials to attempt or solve new loci on the multi-access box (an additional measure of behavioral flexibility), then the **average number of trials to solve** and the **average number of trials to attempt** a new option on the multi-access box will be additional dependent variables. See behavioral flexibility [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md).
3) **Flexibility comprehensive**: This measure is currently being developed and is intended be a more accurate representation of all of the choices an individual made, as well as accounting for the degree of uncertainty exhibited by individuals as preferences change. If this measure more effectively represents flexibility (determined using a modeled dataset and not the actual data), we may decide to solely rely on this measure and not use independent variables 1-3. If this ends up being the case, we will modify the code in the analysis plan below to reflect this change before conducting analyses of the data in this preregistration.
All models will be run once per dependent variable.
**Independent variables**
1) P1: Latency to enter a novel environment inside a familiar environment
2) P1: Time spent in each of the different sections inside a novel environment or the corresponding areas on the floor when the novel environment is not present (familiar environment) as an interaction with the Environment Condition: activity in novel environment vs. activity in familiar environment
3) P1: Time spent per section of a novel environment or in the corresponding areas on othe floor when the novel environment is not present (familiar environment) as an interaction with the Environment Condition: time spent in novel environment vs. time spent in familiar environment
4) P1: Time spent exploring the outside of the novel environment (within 20cm) before entering it
5) P2: Latency to land on the table next to an object (novel, familiar) (that does not contain food) in a familiar environment (that contains maintenance diet away from the object) - OR - latency to touch an object (novel, familiar) (choose the variable with the most data)
6) P3: Number of touches to the functional part of an apparatus per time (multi-access box, novel environment in P1, novel objects in P2 and P4)
7) P3: Number of touches to the non-functional part of an apparatus per time ([multi-access box](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md))
8) P4: Latency to land on the table - OR - Latency to eat the food - OR - Latency to touch a threatening object next to food (choose the variable with the most data)
9) P5: Number of different motor actions used when attempting to solve the [multi-access box](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md)
10) Age (adult: after hatch year, juvenile: hatch year). NOTE: this variable will be removed if only adults are tested (and we are planning to test only adults).
11) ID (random effect because multiple measures per individual)
12) [Condition](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md): control, flexibility manipulation
**ANALYSIS 3** - *H1: P1 alternative 4: correlation between boldness and exploration*
Dependent variable: Boldness: Latency to land on the table - OR - Latency to eat the food - OR - Latency to touch a threatening object next to food (we will choose the variable with the most data)
Independent variables:
1) Time spent exploring the outside of the novel environment (within 20cm) before entering it
2) Latency to land on the table next to an object (novel, familiar) (that does not contain food) in a familiar environment (that contains maintenance diet away from the object) - OR - latency to touch an object (novel, familiar) (choose the variable with the most data)
**ANALYSIS 4** - *H1: P3: does persistence correlate with reversal persistence?*
Dependent variable: The number of incorrect choices in the final reversal before making the first correct choice
Independent variables:
1) Average number of touches to the functional part of an apparatus per time ([multi-access box](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md), novel environment in P1, novel objects in P2 and P4)
2) Condition: control, flexibility manipulation
**ANALYSIS 5** - *H2: P6: captive vs wild*
**Dependent variables**
1) Boldness: In captivity we will measure boldness as the latency to land on the table - OR - Latency to eat the food - OR - Latency to touch a threatening object that is next to food (we will choose the variable with the most data); In the wild the dependent variable will be the latency to come within 2m - OR - Latency to eat the food - OR - Latency to touch a threatening object that is next to food (we will choose the variable with the most data).
2) Persistence: Number of touches to an apparatus per time (multi-access box in the behavioral flexibility [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md), novel environment in P1, objects in P2 and P4)
3) Exploration of novel environment: Latency to enter a novel sub-environment inside a familiar environment
4) Exploration of novel object: Latency to land next to an object (novel, familiar) (that does not contain food) in a familiar environment (that contains maintenance diet away from the object) - OR - latency to touch an object (novel, familiar) (choose the variable with the most data)
*Note: if 3 and 4 are consistent within individuals, and correlate, we will combine these variables into one exploration propensity score.*
**Independent variables**
1) Context: captive or wild
2) Number of times we attempted to assay boldness or exploration but failed due to lack of participation
3) ID (random effect because multiple measures per individual)
### E. ANALYSIS PLAN
We do not plan to **exclude** any data. When **missing data** occur, the existing data for that individual will be included in the analyses for the tests they completed. Analyses will be conducted in R (current version `r getRversion()`; [@rcoreteam]). When there is more than one experimenter within a test, experimenter will be added as a random effect to account for potential differences between experimenters in conducting the tests. If there are no differences between models including or excluding experimenter as a random effect, then we will use the model without this random effect for simplicity.
#### *Unregistered analysis: interobserver reliability of dependent variables*
To determine whether the experimenter coded the dependent variables in a repeatable way, hypothesis-blind video coders were first trained in video coding the dependent variable, and then he coded 26% of the videos in the exploration and boldness experiments. We randomly chose four (Tomatillo, Queso, Mole, and Habanero) of the 19 birds (21%) who participated in these experiments using random.org. Video coders then analyzed all videos from these four birds. The experimenter's data was compared with video coder data using the intra-class correlation coefficient (ICC) to determine the degree of bias in the regression slope (@hutcheon2010random, using the irr package in R: @gamer2012package).
##### Interobserver reliability training
To pass **interobserver reliability (IOR) training**, video coders needed an ICC score of 0.90 or greater to ensure the instructions were clear and that there was a high degree of agreement across coders (see R code comments for details).
**Sierra Planck (discussed with Logan):** Persistence (total number of touches to apparatus) and motor diversity (presence or absence of a behavior from the ethogram). Planck was the first to code videos for these variables so there was not an already established training process or someone to compare her to. Planck and Logan worked together to agree on coding decisions using one video, and then Planck proceeded to code videos independently after that.
**Alexis Breen**
- **Persistence (compared with Logan):** total number of functional touches to apparatus unweighted Cohen's kappa = 1.00 (confidence boundaries=1.00-1.00, n=21 data points)
- **Persistence (compared with Logan):** total number of non-functional touches to apparatus unweighted Cohen's kappa = 0.00 (confidence boundaries=0.00-0.00, n=19 data points). Note: Breen was previously unclear about when to count non-functional touches, however, a discussion eliminated confusion and we proceeded with allowing her to video code independently because the functional touches, which she scored perfectly on, are the more difficult touches to code and thus indicative of her ability to code non-functional touches after clarity on the instructions.
- **Motor diversity (compared with Planck):** presence or absence of a behavior from the ethogram unweighted Cohen's kappa = 0.70 (confidence boundaries=0.39-1.00, n=21 data points). Note: Breen joined the project after Planck and had extensive experience with video coding bird behaviors. Because of this, and because she became Kiepsch's supervisor for exploration, boldness, persistence, and motor diversity, we decided to use Breen as the baseline for persistence and motor diversity and match future coders to her rather than to Plank. Therefore, we moved Breen into the primary video coder position (coding more of the videos than the others). To prepare for Kiepsch's training, Breen clarified the motor diversity ethogram to make it more repeatable. However, we did not require Planck to redo training because she was already so far through the videos. As such, we realize that Planck's data from 21% of the videos may not match Breen's as closely as if Plank was matched to Breen during training.
***Vincent Kiepsch*** (compared with Breen):
- **Exploration** order of the latency-distance categories ICC = 0.96 (confidence boundaries=0.92-1.00, n=141 data points)
- **Boldness** order of the latency-distance categories ICC=1.00 (confidence boundaries=1.00-1.00, n=11 data points). Note that, for exploration and boldness, the ordered categories were aligned based on similar latencies between coders to prevent disagreements near the top of the data sheet from misaligning all subsequent entries.
- **Persistence** number of touches to the apparatus ICC = 0.999 (confidence boundaries=0.996-1.00, n=5 data points).
- **Motor diversity**: the training score for the presence or absence of a behavior from the ethogram required additional training than originally planned, resulting in a final Cohen's kappa = 0.93 (confidence boundaries=0.80-1.00, n=42 data points).
##### Interobserver reliability scores were as follows (4/19 birds; 21% of the videos):
**Vincent Kiepsch** (compared with Breen):
- **Exploration:** closest distance category to apparatus Cohen's unweighted kappa = 0.86 (confidence boundaries=0.71-1.00, n=32 data points)
- **Exploration environment:** first latency to enter tent ICC = 0.997 (confidence boundaries=0.99-0.999, n=10 data points)
- **Boldness:** closest distance to apparatus Cohen's unweighted kappa = 0.86 (confidence boundaries=0.68-1.00, n=24 data points)
**Exploration and boldness in the WILD** (comparison between McCune video coding and transcribing field notes for 20% of the grackles in the wild sample in March 2021 and again on the same data in May 2021):
- Exploration and boldness data collected in the wild were combined because there was not much data for either and because the variables were the same for both assays
- **Exploration and boldness:** closest distance category to apparatus Cohen's unweighted kappa = 1.00 (confidence boundaries=1.00-1.00, n=12 data points)
- **Exploration and boldness:** latency to first landing in a distance category ICC = 0.999 (confidence boundaries=0.994-1.000, n=8 data points)
**Persistence and Motor Diversity** (comparisons between Breen, Kiepsch, and Planck):
- **Persistence:**
- total number of FUNCTIONAL touches to apparatus ICC = 0.77 (confidence boundaries=0.48-0.90, n=18 data points)
- total number of NON-FUNCTIONAL touches to apparatus ICC = 0.68 (confidence boundaries=-0.06-0.95, n=6 data points)
- **Motor diversity:** presence or absence of a behavior from the ethogram unweighted Kappa = 0.77 (confidence boundaries=0.70-0.84, n=380 data points)
These scores indicate that the dependent variables are repeatable to a moderate (persistence and motor diversity) or a high to very high (exploration and boldness) degree given our instructions and training.
```{r ior, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
# Inter/intra-rater reliability using Cohen's kappa when the variable is categorical (scale=1+) or intra-class correlation coefficient when the variable is continuous (Mandrekar 2011 J Thoracic Oncology 6(1):6-7 https://doi.org/10.1097/JTO.0b013e318200f983)
# Intra-class correlation / reliability coefficient / the
# degree of bias in the regression slope (Hutcheon et al.
# 2010. Random measurement error and regression dilution bias
# www.bmj.com/content/340/bmj.c2289). 'The ratio of variation
# in error-free (true) X values to the variation in the
# observed error-prone (observed) values is known as the
# reliability coefficient, attenuation factor, or intra-class
# correlation.'
#Cohen's kappa = Good for nominal data (where distance doesn't mean anything; don't use the weighted Kappa bc it is like the ICC) https://www.rdocumentation.org/packages/psych/versions/1.9.12.31/topics/cohen.kappa
# ICC / Cohen's Kappa must be 0.90 or greater to be considered reliable and pass training
### ICCs for agreement between the 2 coders (live coder and video coder)
#### PASSING interobserver reliability TRAINING so they can become second coders for experiments
#Note: this data counts as second coder data if they have ICC or Kappa > 0.89
#coder 1 = Alexis Breen
#coder 2 = Vincent Kiepsch
library(irr) #ICC package
library(psych) #Cohen's kappa package
##### EXPLORATION
#### IOR TRAINING: Did Vincent pass interobserver reliability TRAINING? YES (5 Nov 2020)
#video file names analyzed: A035P- 2018-12-27 ExploreEnv Time1 Novel T1, A035P- 2019-2-12 ExploreEnv Time2 Novel T1, A035P- 2018-12-28 ExploreObj Time1 Familiar T1
#For each coder, list the order of the distance categories scored (rows aligned to match latencies)
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_expbolmotper_vincent_training.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
de <- d[d$Experiment=="Exploration",]
head(de) #Check to make sure it looks right
de[,3] #coder1
de[,4] #coder2
cohen.kappa(de[,c(3,4)], w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted kappa=0.96, confidence boundaries=0.92-1.00, n=141 data points
#### IOR: On 20% of the videos, what was the interrater reliability between Vincent and Alexis?
# Exploration OBJECT and ENVIRONMENT: Closest distance to apparatus
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_exploration_vincent.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
d[,4] #Coder1closestdistancetoapparatus
d[,6] #Coder2closestdistancetoapparatus
cohen.kappa(d[,c(4,6)], w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
# unweighted kappa=0.86, confidence boundaries: 0.71-1.00, n=32 data points
# Exploration ENVIRONMENT: First latency to enter tent
d[,3] #Coder1timeentertent
d[,5] #Coder2timeentertent
# Need to change column 3 & 5 to numeric (b.c currently categorical)
d$Coder1timeentertent<-as.numeric(d$Coder1timeentertent)
d$Coder2timeentertent<-as.numeric(d$Coder2timeentertent)
icc(d[,c(3,5)], model="oneway", type="consistency", unit="single", conf.level=0.95)
# ICC = 0.997 (95% Confidence Interval: 0.99-0.999, n=32 data points (but only 10 with data in them))
##### BOLDNESS
#### TRAINING: Did Vincent pass interobserver reliability TRAINING? YES
#video file names analyzed: A022-N 2019-01-15 Boldness Time2 C2 T2, A031-Y 2018-12-24 Boldness C2 T2
#For each coder, list the order of the distance category per video
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_expbolmotper_vincent_training.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
db <- d[d$Experiment=="Boldness",]
head(db) #Check to make sure it looks right
db[,3]
db[,4]
cohen.kappa(db[,c(3,4)], w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted kappa=1.00, confidence boundaries=1.00-1.00, n=11 data points
#### IOR: On 20% of the videos, what was the interrater reliability between Vincent and Alexis?
##The data for IOR_explorationboldness_vincent.csv comes from GrackleVideoCoding, tab Exp7&8_IOR_AB_VK
# Closest distance to apparatus
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_exploration_vincent.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
db <- d[d$Experiment=="Boldness",]
head(db) #Check to make sure it looks right
db[,4] #coder1
db[,6] #coder2
db$Coder1closestdistancetoapparatus<-as.factor(db$Coder1closestdistancetoapparatus)
db$Coder2closestdistancetoapparatus<-as.factor(db$Coder2closestdistancetoapparatus)
#load tidyverse for data wrangling
library(tidyverse)
#clean col 4 & 6 b/c if not than counts n/a as a factor level in IOR calc
dbclean<-db %>%
filter(Coder1closestdistancetoapparatus != "n/a") %>% #throw out n/a category
filter(Coder2closestdistancetoapparatus != "n/a") %>% #throw out n/a category
droplevels() #drop unused level n/a from factor 4 & 6
cohen.kappa(dbclean[,c(4,6)], w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
# unweighted kappa=0.86, confidence boundaries: 0.68-1.00, n=24 data points
#### WILD BOLDNESS & EXPLORATION ####
#Coded by the experimenter 1) live in the field using video and field notes, 2) more than a year later from video and field notes, and 3) 20% of individuals re-coded from video and field notes more than 2 months later for IOR calculation.
fd <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_expbolfieldKM.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(fd) #Check to make sure it looks right
## Exploration AND boldness (not much data so they are combined): latency to land in a distance category
fd[,5] #Coder1latencyland
fd[,7] #Coder2latencyland
# Need to change columns to numeric (b.c currently categorical)
fd$Coder1latencyland<-as.numeric(fd$Coder1latencyland)
fd$Coder1latencyland<-as.numeric(fd$Coder2latencyland)
icc(fd[,c(5,7)], model="oneway", type="consistency", unit="single", conf.level=0.95)
# ICC = 0.999 (95% Confidence Interval: 0.994-1.000, n=8 data points (but only 10 with data in them))
## Exploration AND boldness (not much data so they are combined): closest distance category
fd[,6]
fd[,8]
cohen.kappa(fd[,c(6,8)], w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted kappa=1.00, confidence boundaries=1.00-1.00, n=12 data points
##### MOTOR DIVERSITY
#coder 1 = Corina Logan
#coder 2 = Sierra Planck
#coder 3 = Alexis Breen
#coder 4 = Vincent Kiepsch
#### IOR TRAINING for motor diversity and persistence at the same time: Did Sierra calibrate to Corina? Yes
#Corina and Sierra both coded A035P- 2018-12-12 Multiaccess box S1 T1 to calibrate
# motor diversity:
#Corina=8 motor actions (data sheet: data_IORcorina at Google Drive): PullVert, PeckHorz, Grab, Pull, PeckHorz, GrabUpDown, PullHorz, PeckVert
#Sierra=6 motor actions (GrackleVideoCoding google sheet: tab Exp_mot row 2):
#Corina decided Sierra passed because we went back through the video together. NOTE: Sierra has different codings here and Corina agrees with her - we worked out some kinks after Corina coded this data and we improved the categories. The only thing Sierra didn't have that Corina did was Vertical peck 3:04.
#Decision: Sierra moves forward with video coding bc she is coding 100% and we will match the IOR to her
#Sierra then recoded the same video=9 motor actions (GrackleVideoCoding google sheet: tab Exp_mot row 2): GrabVert, PullVert, PeckHorz, PeckVert, GrabHorz, PeckUpDown, GrabUpDown, PullHorz, PushHorz
# Persistence: Total touches to apparatus: Corina=40, Sierra=40 - perfect agreement
# Matching Alexis to Sierra
#Alexis coded A035P- 2018-12-12 Multiaccess box S1 T1 = 9 motor actions (data sheet: GrackleVideoCoding google sheet: tab IOR_mot row 2). Alexis listed 3 actions (out of 9 total) that Sierra did not list
#Alexis n=9: GrabVert, PullVert, PeckVert (not in Sierra's), PeckHorz, GrabHorz, PeckUpDown, GrabUpDown (not in Sierra's), PullHorz (not in Sierra's), PushHorz
#01:12, 01:12, 01:34, 01:37, 01:45, 02:24, 02:25, 02:39, 02:41
#Compare Alexis and Sierra: 21 motor actions in the ethogram, score presence/absence of each
ab <- c(1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0)
sp <- c(1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
df <- data.frame(ab,sp)
cohen.kappa(df, w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted kappa=0.70, 0.39-1.00, n=21 data points
# Persistence: Total touches to apparatus: Alexis=32. All of us (Alexis, Corina, Sierra) counted n=21 functional touches, the difference was in the nonfunctional touches
#Functional touches: compare CL and AB
cl <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
ab2 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
df2 <- data.frame(cl,ab2)
cohen.kappa(df2, w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted kappa=1.00, 1.00-1.00, n=21 data points
#NON-functional touches: compare CL and AB
cl2 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
ab3 <- c(1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0)
df3 <- data.frame(cl2,ab3)
cohen.kappa(df3, w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted kappa=0.00, 0.00-0.00, n=19 data points
#Decided to have Alexis code all videos for bird 35 and match all future coders to her because she made the motor diversity ethogram more detailed and with diagrams, and has an extensive background in coding these kinds of behaviors. So move Sierra to the 20% person for MAB plastic and MAB log. Vincent is the 20% person for Exp and Bol.
#From now (6 Oct 2020) on, for somone to pass interobserver reliability training, they must get Cohen's kappa 0.9 or above for MotorActionsUsed column motor action order, and ICC 0.9 or above for TotalNumberFunctionalTouches column time stamps
#### IOR TRAINING: Did Vincent pass interobserver reliability TRAINING for motor diversity on exploration & boldness videos? NO (4 Nov 2020), then did more training and did pass
#video file names analyzed 5 Nov 2020: A035P- 2018-12-27 ExploreEnv Time1 Novel T1
#video file names analyzed 6 Nov 2020: A035P- 2018-12-27 ExploreEnv Time1 Novel T1, A035P- 2019-2-12 ExploreEnv Time2 Novel T1, A035P- 2018-12-28 ExploreObj Time1 Familiar T1, A022-N 2019-01-15 Boldness Time2 C2 T2, A031-Y 2018-12-24 Boldness C2 T2, A022-N 2018-12-10 ExploreEnv Novel T1
#video file names analyzed 11 Nov 2020 (analysis 1): A035P- 2018-12-27 ExploreEnv Time1 Novel T1, A035P- 2019-2-12 ExploreEnv Time2 Novel T1, A035P- 2018-12-28 ExploreObj Time1 Familiar T1, A022-N 2019-01-15 Boldness Time2 C2 T2, A031-Y 2018-12-24 Boldness C2 T2, A022-N 2018-12-10 ExploreEnv Novel T1, A022-N 2019-01-11 ExploreEnv Time2 Novel T1
#video file names analyzed 11 Nov 2020 (analysis 2 without the first 2 videos=38 data points): A035P- 2018-12-28 ExploreObj Time1 Familiar T1, A022-N 2019-01-15 Boldness Time2 C2 T2, A031-Y 2018-12-24 Boldness C2 T2, A022-N 2018-12-10 ExploreEnv Novel T1, A022-N 2019-01-11 ExploreEnv Time2 Novel T1
#video file names analyzed 11 Nov 2020 (analysis 3 only bird 22 explore environment): A022-N 2018-12-10 ExploreEnv Novel T1, A022-N 2019-01-11 ExploreEnv Time2 Novel T1
#video file names analyzed 11 Nov 2020: A035P- 2018-12-27 ExploreEnv Time1 Novel T1, A035P- 2019-2-12 ExploreEnv Time2 Novel T1, A035P- 2018-12-28 ExploreObj Time1 Familiar T1, A022-N 2019-01-15 Boldness Time2 C2 T2, A031-Y 2018-12-24 Boldness C2 T2, A022-N 2018-12-10 ExploreEnv Novel T1, A022-N 2019-01-11 ExploreEnv Time2 Novel T1, A073OL 2019-11-28 ExploreEnv Time2 Novel T1 (pt. 1)
#video file names analyzed 16 Nov 2020: A035P- 2018-12-27 ExploreEnv Time1 Novel T1, A035P- 2019-2-12 ExploreEnv Time2 Novel T1, A035P- 2018-12-28 ExploreObj Time1 Familiar T1, A022-N 2019-01-15 Boldness Time2 C2 T2, A031-Y 2018-12-24 Boldness C2 T2, A022-N 2018-12-10 ExploreEnv Novel T1, A022-N 2019-01-11 ExploreEnv Time2 Novel T1, A073OL 2019-11-28 ExploreEnv Time2 Novel T1 (pt. 1), A073OL 2019-10-25 ExploreEnv Time1 Novel T1
#video file names analyzed 16 Nov 2020 (only those videos after the intensive training session): A073OL 2019-11-28 ExploreEnv Time2 Novel T1 (pt. 1), A073OL 2019-10-25 ExploreEnv Time1 Novel T1 = PASSED
# motor diversity: presence/absence of a coder coding a particular behavior from the ethogram
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_expbolmotper_vincent_training.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
dm <- d[d$Experiment=="Motor",]
head(dm) #Check to make sure it looks right
dm[,3]
dm[,4]
cohen.kappa(dm[,c(3,4)], w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted = 0.74, 0.61-0.87, n=105 data points. Did not pass on 5 Nov 2020
#unweighted = 0.74, 0.62-0.86, n=126 data points. Did not pass on 6 Nov 2020
#unweighted = 0.75, 0.64-0.87, n=147 data points. Did not pass on 11 Nov 2020 (1). There were too many mistakes in the first two videos for us to be able to evaluate whether he now passes training. Therefore, I am removing the first two videos from the analysis and recalculating...
#unweighted = 0.74, 0.55-0.92, n=109 data points. Did not pass on 11 Nov 2020 (2). Try only analyzing the most recent videos coded on bird 22 (explore env), which occurred after more training
#unweighted = 0.69, 0.41-0.99, n=42 data points. Did not pass on 11 Nov 2020 (3). So even his recent training is not improving his scoring ability. It's better to include all of the videos coded so far.
#unweighted = 0.8, 0.71-0.89, n=168 data points. Did not pass on 13 Nov 2020
#unweighted = 0.82, 0.73-0.90, n=189 data points. Did not pass on 16 Nov 2020
#unweighted = 0.93, 0.80-1.00, n=42 data points. PASSED on 16 Nov 2020 (need to replace the other training videos to make sure he codes 20% of all videos)
#### On 20% of the videos, what was the interrater reliability between Vincent and Alexis?
#total number of motor actions per bird across all experiments
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_motor_BreenPlanckKiepsch.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
d[,6]
d[,7]
cohen.kappa(d[,c(6,7)], w=NULL,n.obs=NULL,alpha=.05,levels=NULL)
#unweighted kappa 0.77 (0.70-0.84), n`380 data points
##### PERSISTENCE
#coder 1=Alexis Breen
#coder 2=Vincent Kiepsch
#### IOR TRAINING: Did Vincent pass interobserver reliability TRAINING for persistence on exploration & boldness videos? YES (4 Nov 2020)
#video file names analyzed: A035P- 2018-12-27 ExploreEnv Time1 Novel T1, A035P- 2019-2-12 ExploreEnv Time2 Novel T1, A035P- 2018-12-28 ExploreObj Time1 Familiar T1, A022-N 2019-01-15 Boldness Time2 C2 T2, A031-Y 2018-12-24 Boldness C2 T2
#list the values in the TotalNumberFunctionalTouches for each row (a row is a video or a trial depending on the experiment) (note: for exploration and boldness, there are only functional touches and no nonfunctional touches)
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_expbolmotper_vincent_training.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
dp <- d[d$Experiment=="Persistence",]
head(dp) #Check to make sure it looks right
dp[,3] <- as.numeric(dp[,3]) #coder1
dp[,4] <- as.numeric(dp[,4]) #coder2
icc(dp[,c(3,4)], model = "oneway", type = "consistency", unit = "single", conf.level = 0.95)
#ICC=0.999, 95%-Confidence Interval 0.996 < ICC < 1, n=5 data points
#### IOR: On 20% of the videos, what was the interrater reliability between Vincent and Alexis?
#total number of touches to an apparatus per bird across all experiments
#need to change the file path to the correect tab at the google sheet - the one here is for training data
d <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/IOR_persistence_BreenPlanckKiepsch.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(d) #Check to make sure it looks right
d[,3] <- as.numeric(d[,3]) #coder1
d[,4] <- as.numeric(d[,4]) #coder2
d[,5] <- as.numeric(d[,5]) #coder1
d[,6] <- as.numeric(d[,6]) #coder2
icc(d[,c(3,4)], model = "oneway", type = "consistency", unit = "single", conf.level = 0.95)
#functional touches: ICC=0.77, 95%-Confidence Interval 0.484 < ICC < 0.904, n=18 data points, F-Test, H0: r0 = 0 ; H1: r0 > 0 F(17,18) = 7.52 , p = 4.66e-05
icc(d[,c(5,6)], model = "oneway", type = "consistency", unit = "single", conf.level = 0.95)
#nonfunctional touches: ICC = 0.682, 95%-Confidence Interval -0.062 < ICC < 0.947, n=6 data points, F-Test, H0: r0 = 0 ; H1: r0 > 0 F(5,6) = 5.29 , p = 0.0332
```
#### *Ability to detect actual effects*
To begin to understand what kinds of effect sizes we will be able to detect given our sample size limitations and our interest in decreasing noise by attempting to measure it, which increases the number of explanatory variables, we used G&ast;Power (v.3.1, @faul2007g, @faul2009statistical) to conduct power analyses based on confidence intervals. G&ast;Power uses pre-set drop down menus and we chose the options that were as close to our analysis methods as possible (listed in each analysis below). Note that there were no explicit options for GLMs (though the chosen test in G&ast;Power appears to align with GLMs) or GLMMs or for the inclusion of the number of trials per bird (which are generally large in our investigation), thus the power analyses are only an approximation of the kinds of effect sizes we can detect. We realize that these power analyses are not fully aligned with our study design and that these kinds of analyses are not appropriate for Bayesian statistics (e.g., our MCMCglmm below), however we are unaware of better options at this time. Additionally, it is difficult to run power analyses because it is unclear what kinds of effect sizes we should expect due to the lack of data on this species for these experiments.
To address the power analysis issues, we will run simulations on our Arizona data set before conducting any analyses in this preregistration. We will first run null models (i.e., dependent variable ~ 1 + random effects), which will allow us to determine what a weak versus a strong effect is for each model. Then we will run simulations based on the null model to explore the boundaries of influences (e.g., sample size) on our ability to detect effects of interest of varying strengths. If simulation results indicate that our Arizona sample size is not larger than the lower boundary, we will continue these experiments at the next field site until we meet the minimum suggested sample size.
#### *Data checking*
The data will be checked for overdispersion, underdispersion, zero-inflation, and heteroscedasticity with the DHARMa R package [@Hartig2019dharma] following methods by [Hartig](https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html). Note: DHARMa doesn't support MCMCglmm, therefore we will use the closest supported model: glmer from the R package lme4 [@lme4].
#### *Repeatability of exploration, boldness and persistence*
**Analysis:** We will obtain repeatability estimates that account for the observed and latent scales, and then compare them with the raw repeatability estimate from the null model. The repeatability estimate indicates how much of the total variance, after accounting for fixed and random effects, is explained by individual differences (ID). We will run this GLMM using the MCMCglmm function in the MCMCglmm package ([@hadfieldMCMCglmmpackage]) with a Poisson distribution and log link using 13,000 iterations with a thinning interval of 10, a burnin of 3,000, and minimal priors (V=1, nu=0) [@hadfield2014coursenotes]. We will ensure the GLMM shows acceptable convergence (i.e., lag time autocorrelation values <0.01; [@hadfieldMCMCglmmpackage]), and adjust parameters if necessary.
Note Feb 2021: a Gaussian distribution was used instead of a Poisson for exploration and boldness latencies because they are continuous variables.
Note: The power analysis is the same as for P3 (below) because there are the same number of explanatory variables (fixed effects).
```{r repeatableBold, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Boldness
bol <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/g_exploration_data_boldness.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(bol) #Check to make sure it looks right
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(LatencyToFeed ~ Condition + (1|ID), family=poisson, data=bol), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#REPEATABILITY
#GLMM
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
bold <- MCMCglmm(LatencyToFeed ~ Condition, random=~ID, family="poisson", data=bol, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(bold)
#autocorr(bold$Sol) #Did fixed effects converge?
#autocorr(bold$VCV) #Did random effects converge?
#Advice from Maxime Dahirel in the review for flexmanip on how to calculate repeatability from MCMCglmm (https://ecology.peercommunityin.org/public/rec?id=17&reviews=True). Dahirel says: “7- In P3a, you use a linear mixed model (lmer) to estimate repeatability when a generalized linear mixed model (glmer) should be used, as the dependent variable is a count. This would lead to potentially wrong estimates of repeatability, first because the individual variance will probably be wrongly estimated, and then because the residual/overdispersion variance is not the same for a Gaussian vs a Poisson model (Nakagawa & Schielzeth, 2010). If you want to extract correct repeatabilities from a lme4 (G)LMM (and their 95% CI) on both the response and latent scales, you can use the rptR package (Stoffel, Nakagawa, & Schielzeth, 2017). However, since you plan to fit your GLMM using MCMCglmm you actually don’t need that at all, since the latent scale adjusted repeatability and its credible interval can simply be obtained by mod$VCV[,ID]/(mod$VCV[,ID]+mod$VCV[,units]+mod$VCV[,any other included random effect]). For the raw repeatability, simply fit a model with no covariates. For the repeatability on the response scale, see (P. de Villemereuil, Morrissey, Nakagawa, & Schielzeth, 2018; Pierre de Villemereuil, Schielzeth, Nakagawa, & Morrissey, 2016) and the QGglmm R package.”
#Therefore, because MCMCglmm can easily calculate repeatability on the latent scale, there is no need to transform this back to the observed scale.
repeata <- bold$VCV[,"ID"]/(bold$VCV[,"ID"]+bold$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
mean(repeata)
var(repeata)
posterior.mode(repeata)
HPDinterval(repeata, 0.95)
# Repeatability on the data/observed scale (accounting for fixed effects)
#code from Supplementary Material S2 from Villemereuil et al. 2018 J Evol Biol
vf <- sapply(1:nrow(bold[["Sol"]]), function(i) {
var(predict(bold, it=i))
}) #estimates for each iteration of the MCMC
repeataF <- (vf+bold$VCV[,"ID"])/(vf+bold$VCV[,"ID"]+bold$VCV[,"units"]) #latent scale adjusted + data scale
posterior.mode(repeataF)
HPDinterval(repeataF, 0.95)
# Now compare with the raw repeatability: null model
boldraw <- MCMCglmm(LatencyToFeed ~ 1, random=~ID, family="poisson", data=bol, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(boldraw)
repeataraw <- boldraw$VCV[,"ID"]/(boldraw$VCV[,"ID"]+boldraw$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
posterior.mode(repeata)
HPDinterval(repeata, 0.95)
```
```{r repeatablePer, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Persistence
per <- read.csv ("/Users/corina/GTGR/data/data_persist.csv", header=T, sep=",", stringsAsFactors=F)
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(NoTouches ~ Test * Condition + (1|ID), family=poisson, data=per), n=250)
simulationOutput$scaledResiduals
testDispersion(simulationOutput)
testZeroInflation(simulationOutput)
testUniformity(simulationOutput)
plot(simulationOutput)
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #can't get this code to work yet
#REPEATABILITY
#GLMM
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
pers <- MCMCglmm(NoTouches ~ Test*Condition, random=~ID, family="poisson", data=per, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(pers)
#autocorr(pers$Sol) #Did fixed effects converge?
#autocorr(pers$VCV) #Did random effects converge?
#In MCMCglmm, the latent scale adjusted repeatability and its credible interval can simply be obtained by: mod$VCV[,ID]/(mod$VCV[,ID]+mod$VCV[,units]) - advice from Maxime Dahirel
repeata <- pers$VCV[,"ID"]/(pers$VCV[,"ID"]+pers$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
mean(repeata)
var(repeata)
posterior.mode(repeata)
HPDinterval(repeata, 0.95)
# Repeatability on the data/observed scale (accounting for fixed effects)
#code from Supplementary Material S2 from Villemereuil et al. 2018 J Evol Biol
vf <- sapply(1:nrow(pers[["Sol"]]), function(i) {
var(predict(pers, it=i))
}) #estimates for each iteration of the MCMC
repeataF <- (vf+pers$VCV[,"ID"])/(vf+pers$VCV[,"ID"]+pers$VCV[,"units"]) #latent scale adjusted + data scale
posterior.mode(repeataF)
HPDinterval(repeataF, 0.95)
# Now compare with the raw repeatability: null model
persraw <- MCMCglmm(NoTouches ~ 1, random=~ID, family="poisson", data=per, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(persraw)
repeataraw <- persraw$VCV[,"ID"]/(persraw$VCV[,"ID"]+persraw$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
posterior.mode(repeata)
HPDinterval(repeata, 0.95)
```
```{r repeatableExpEnv, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Exploration of novel environment
ee <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/g_exploration_data_exploration.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(ee) #Check to make sure it looks right
#want only Test=Environment (not Object)
ee <- ee[ee$Test=="Environment",]
#LatencyFirstLand has the most data so will use this variable
#should I change LatencyFirstLand NAs to a number higher than the assay session length so they show up as data points? yes. Set latency for these at 3400 (the highest assay duration was 3340s)
ee$LatencyFirstLand <- as.numeric(ee$LatencyFirstLand)
ee$BirdID <- as.factor(ee$BirdID)
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(log(LatencyFirstLand) ~ Condition + (1|BirdID), family=gaussian, data=ee), n=250) #250 simulations, but if want higher precision change n>1000
plot(simulationOutput$scaledResiduals) #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
#looks randomly spread
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
#p=0.832 so not over or under dispersed
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p>0.05
#p=1 so not zero inflated
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05.
# p = 0.035, data not normally distributed.
# Log-transformed Latency variable, and now p-value = 0.23
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ee$Test, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#REPEATABILITY
#GLMM - 8 Feb 2021: changed family from poisson to gaussian because latency is gaussian, have to log-transform Latency so residuals are normally distributed
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
expl <- MCMCglmm(log(LatencyFirstLand) ~ Condition, random=~BirdID, family="gaussian", data=ee, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(expl)
autocorr(expl$Sol) #Did fixed effects converge? Values after Lag 0 should be < 0.1. Most values=Yes
autocorr(expl$VCV) #Did random effects converge? Most values=Yes
#Advice from one of our reviewers Maxime Dahirel: In MCMCglmm, the latent scale adjusted repeatability and its credible interval can simply be obtained by: mod$VCV[,ID]/(mod$VCV[,ID]+mod$VCV[,units])
repeata <- expl$VCV[,"BirdID"]/(expl$VCV[,"BirdID"]+expl$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
mean(repeata) #0.37
var(repeata) #0.02
posterior.mode(repeata) #0.33
HPDinterval(repeata, 0.95) #0.10-0.63
# Repeatability on the data/observed scale (accounting for fixed effects)
#code from Supplementary Material S2 from Villemereuil et al. 2018 J Evol Biol
vf <- sapply(1:nrow(expl[["Sol"]]), function(i) {
var(predict(expl, it=i))
}) #estimates for each iteration of the MCMC
repeataF <- (vf+expl$VCV[,"BirdID"])/(vf+expl$VCV[,"BirdID"]+expl$VCV[,"units"]) #latent scale adjusted + data scale
posterior.mode(repeataF) #0.35
HPDinterval(repeataF, 0.95) #0.11-0.63
# Now compare with the raw repeatability: null model
explraw <- MCMCglmm(log(LatencyFirstLand) ~ 1, random=~BirdID, family="gaussian", data=ee, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(explraw)
repeataraw <- explraw$VCV[,"BirdID"]/(explraw$VCV[,"BirdID"]+explraw$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
posterior.mode(repeata) #0.33
HPDinterval(repeata, 0.95) #0.09-0.63
#RESULT: Exp Env = repeatable
## Unregistered analysis - compare R values with rptR function and LRT test to make sure the interpretation of MCMCglmm results is correct - output matches that of MCMCglmm
library(rptR)
r1 = rpt(log(LatencyFirstLand) ~ Condition + (1|BirdID), grname = "BirdID", data=ee, datatype="Gaussian", nboot=1000, npermut=100)
summary(r1)
# R = 0.39
# CI = 0.12-0.60
# p < 0.05
m1 = glmer(log(LatencyFirstLand) ~ Condition + (1|BirdID), data = ee, family = gaussian)
m2 = glm(log(LatencyFirstLand) ~ Condition, data = ee, family = gaussian)
anova(m1,m2) #random effect improves fit, accounts for a significant proportion of the variance.
library(ggplot2)
ggplot(ee[which(ee$Condition == "Novel"),], aes(x = as.factor(Time), y = LatencyFirstLand, colour = BirdID)) +
geom_point() +
facet_wrap( ~ BirdID)
```
```{r repeatableExpObj, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Exploration of novel object
eo <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/g_exploration_data_exploration.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
head(eo) #Check to make sure it looks right
#want only Test=Object (not Environment)
eo <- eo[eo$Test=="Object",]
#LatencyFirstLand has the most data so will use this variable (=3400s when a bird never came to the ground; the highest assay duration was 3340s)
eo$LatencyFirstLand <- as.numeric(eo$LatencyFirstLand)
eo$BirdID <- as.factor(eo$BirdID)
eo$Condition <- as.factor(eo$Condition)
# DATA CHECKING - changed family to gaussian (from poisson) bc latency is gaussian
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(LatencyFirstLand ~ Condition + (1|BirdID), family=gaussian, data=eo, na.action = na.omit ), n=250) #250 simulations, but if want higher precision change n>1000
plot(simulationOutput$scaledResiduals) #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
#looks randomly scattered
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
#p=0.984 so not over or underdispersed
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p>0.05; p = 1
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05
#p=1 so not heteroscadistic
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(eo$Condition, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#REPEATABILITY
#GLMM
#eo2 <- subset(eo,!(is.na(eo["Condition"]))) *KM does not see NA values in the condition column of eo
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
explo <- MCMCglmm(LatencyFirstLand ~ Condition, random=~BirdID, family="gaussian", data=eo, verbose=F, prior=prior, nitt=630000, thin=1000, burnin=23000)
summary(explo) #slower to approach the novel object compared with the familiar object
autocorr(explo$Sol) #Did fixed effects converge? Yes=all values after lag 0 are <0.1 indicating they converged
autocorr(explo$VCV) #Did random effects converge? Yes=two values >0.1, increased nitt (to 630,000) and thin (to 1,000), then all values <0.1
#In MCMCglmm, the latent scale adjusted repeatability and its credible interval can simply be obtained by: mod$VCV[,ID]/(mod$VCV[,ID]+mod$VCV[,units]) - advice from Maxime Dahirel
repeata <- explo$VCV[,"BirdID"]/(explo$VCV[,"BirdID"]+explo$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
mean(repeata) #0.01
var(repeata) #0.002
posterior.mode(repeata) #-1.03e-06
HPDinterval(repeata, 0.95) #7.91e-23 to 0.08
# Repeatability on the data/observed scale (accounting for fixed effects)
#code from Supplementary Material S2 from Villemereuil et al. 2018 J Evol Biol
vf <- sapply(1:nrow(explo[["Sol"]]), function(i) {
var(predict(explo, it=i))
}) #estimates for each iteration of the MCMC
repeataF <- (vf+explo$VCV[,"BirdID"])/(vf+explo$VCV[,"BirdID"]+explo$VCV[,"units"]) #latent scale adjusted + data scale
posterior.mode(repeataF) #0.29
HPDinterval(repeataF, 0.95) #0.11-0.43
# Now compare with the raw repeatability: null model
exploraw <- MCMCglmm(LatencyFirstLand ~ 1, random=~BirdID, family="gaussian", data=eo, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(explraw)
repeataraw <- exploraw$VCV[,"BirdID"]/(exploraw$VCV[,"BirdID"]+exploraw$VCV[,"units"]) #latent scale adjusted repeatability and its credible interval
posterior.mode(repeata) #-1.03e-06
HPDinterval(repeata, 0.95) #7.91e-23 to 0.08
#RESULT=Env Obj is not repeatable
## Unregistered analysis - compare R values with rptR function and LRT test to make sure the interpretation of MCMCglmm results is correct
library(rptR)
r2 = rpt(LatencyFirstLand ~ Condition + (1|BirdID), grname = "BirdID", data=eo, datatype="Gaussian", nboot=1000, npermut=100)
summary(r2)
# R = 0.11
# CI = 0-0.34
# p > 0.05
m3 = glmer(LatencyFirstLand ~ Condition + (1|BirdID), data = eo[-which(is.na(eo$LatencyFirstLand)),], family = gaussian)
m4 = glm(LatencyFirstLand ~ Condition, data = eo[-which(is.na(eo$LatencyFirstLand)),], family = gaussian)
anova(m3,m4) #random effect does not improve the fit of the model, does not account for a significant portion of the variance.
library(ggplot2)
ggplot(eo[which(eo$Condition == "Novel"),], aes(x = as.factor(Time), y = LatencyFirstLand, colour = BirdID)) +
geom_point() +
facet_wrap( ~ BirdID)
```
#### *H1: P1-P5: correlation of flexibility with exploration of new environments and objects, boldness, persistence, and motor diversity*
**Analysis:** If behavior is not repeatable across assays at Time 1 and Time 2 (six weeks apart, both assays occur after the flexibility manipulation takes place) for exploration, boldness, persistence, or motor diversity (see analysis for P6), we will not include these variables in analyses involving flexibility. If behavior is repeatable within individuals, we will examine the relationship between flexibility and these variables as follows. Note that the two exploration measures (novel environment and novel object) will be combined into one variable if they correlate and are both repeatable within individuals.
Because the independent variables could influence each other, we will analyze them in a single model: Generalized Linear Mixed Model (GLMM; MCMCglmm function, MCMCglmm package; [@hadfieldMCMCglmmpackage]) with a Poisson distribution and log link using 13,000 iterations with a thinning interval of 10, a burnin of 3,000, and minimal priors (V=1, nu=0) [@hadfield2014coursenotes]. We will ensure the GLMM shows acceptable convergence (i.e., lag time autocorrelation values <0.01; [@hadfieldMCMCglmmpackage]), and adjust parameters if necessary. We will determine whether an independent variable had an effect or not using the Estimate in the full model.
To roughly estimate our ability to detect actual effects (because these power analyses are designed for frequentist statistics, not Bayesian statistics), we ran a power analysis in G&ast;Power with the following settings: test family=F tests, statistical test=linear multiple regression: Fixed model (R^2 deviation from zero), type of power analysis=a priori, alpha error probability=0.05. We reduced the power to 0.70 and increased the effect size until the total sample size in the output matched our projected sample size (n=32). The number of predictor variables was restricted to only the fixed effects because this test was not designed for mixed models. The protocol of the power analysis is here:
*Input:*
Effect size f² = 0,62
α err prob = 0,05
Power (1-β err prob - note: β=probability of making a Type II error) = 0,7
Number of predictors = 10
*Output:*
Noncentrality parameter λ = 19,8400000
Critical F = 2,3209534
Numerator df = 10
Denominator df = 21
Total sample size = 32
Actual power = 0,7027626
This means that, with our sample size of 32, we have a 70% chance of detecting a large effect (approximated at f^2=0.35 by @cohen1988statistical).
```{r explore, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
explore <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/g_exploration_data_aviaryexplorationmotpersist.csv", header = T, sep = ",", stringsAsFactors = F)
# DATA CHECKING for 1st GLMM
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(TrialsToReverseLast ~ Condition + TimeOutsideNovelEnv + LatencyExpEnv + AverageTimePerSection*EnvCondition + TotalNumberSections*EnvCondition + LatencyTableExpObject + MultiaccessTouchesPerTime + LatencyBoldness + NoMotorActions + (1|Batch), family=poisson, data=explore), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
# ANALYSIS
#Take the average of Time 1 and Time 2 for each variable (exploration environment, exploration object, boldness, motor diversity, persistence)
TimeOutsideNovelEnvT1 <- TimeOutsideNovelEnv[TimeOutsideNovelEnv$Time==1,]
TimeOutsideNovelEnvT2 <- TimeOutsideNovelEnv[TimeOutsideNovelEnv$Time==2,]
TimeOutsideNovelEnv <- (TimeOutsideNovelEnvT1 + TimeOutsideNovelEnvT2)/2
LatencyExpEnvT1 <- LatencyExpEnv[LatencyExpEnv$Time==1,]
LatencyExpEnvT2 <- LatencyExpEnv[LatencyExpEnv$Time==2,]
LatencyExpEnv <- (LatencyExpEnvT1 + LatencyExpEnvT2)/2
AverageTimePerSectionNovelEnvT1 <- AverageTimePerSectionNovelEnv[AverageTimePerSectionNovelEnv$Time==1,]
AverageTimePerSectionNovelEnvT2 <- AverageTimePerSectionNovelEnv[AverageTimePerSectionNovelEnv$Time==2,]
AverageTimePerSectionNovelEnv <- (AverageTimePerSectionNovelEnvT1 + AverageTimePerSectionNovelEnvT2)/2
TotalNumberSectionsNovelEnvT1 <- TotalNumberSectionsNovelEnv[TotalNumberSectionsNovelEnv$Time==1,]
TotalNumberSectionsNovelEnvT2 <- TotalNumberSectionsNovelEnv[TotalNumberSectionsNovelEnv$Time==2,]
TotalNumberSectionsNovelEnv <- (TotalNumberSectionsNovelEnvT1 + TotalNumberSectionsNovelEnvT2)/2
LatencyTableExpObjectT1 <- LatencyTableExpObject[LatencyTableExpObject$Time==1,]
LatencyTableExpObjectT2 <- LatencyTableExpObject[LatencyTableExpObject$Time==2,]
LatencyTableExpObject <- (LatencyTableExpObjectT1 + LatencyTableExpObjectT2)/2
MultiaccessTouchesPerTimeT1 <- MultiaccessTouchesPerTime[MultiaccessTouchesPerTime$Time==1,]
MultiaccessTouchesPerTimeT2 <- MultiaccessTouchesPerTime[MultiaccessTouchesPerTime$Time==2,]
MultiaccessTouchesPerTime <- (MultiaccessTouchesPerTimeT1 + MultiaccessTouchesPerTimeT2)/2
LatencyBoldnessT1 <- LatencyBoldness[LatencyBoldness$Time==1,]
LatencyBoldnessT2 <- LatencyBoldness[LatencyBoldness$Time==2,]
LatencyBoldness <- (LatencyBoldnessT1 + LatencyBoldnessT2)/2
NoMotorActionsT1 <- NoMotorActions[NoMotorActions$Time==1,]
NoMotorActionsT1 <- NoMotorActions[NoMotorActions$Time==1,]
NoMotorActions <- (NoMotorActionsT1 + NoMotorActionsT2)/2
#GLMM - dependent variable = number of trials to reverse in the last reversal
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0),R2=list(V=1,nu=0),R3=list(V=1,nu=0),R4=list(V=1,nu=0),R5=list(V=1,nu=0),R6=list(V=1,nu=0),R7=list(V=1,nu=0),R8=list(V=1,nu=0),R9=list(V=1,nu=0),R10=list(V=1,nu=0),R11=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0),G2=list(V=1,nu=0)))
expl1 <- MCMCglmm(TrialsToReverseLast ~ Condition + TimeOutsideNovelEnv + LatencyExpEnv + AverageTimePerSection*EnvCondition + TotalNumberSections*EnvCondition + LatencyTableExpObject + MultiaccessTouchesPerTime + LatencyBoldness + NoMotorActions, random=~Batch, family="poisson", data=explore, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(expl1)
autocorr(expl1$Sol) #Did fixed effects converge?
autocorr(expl1$VCV) #Did random effects converge?
```
#### *H1: P1 alternative 4: correlations between exploration and boldness measures*
**Analysis:** Generalized Linear Model (GLM; glm function, stats package) with a Poisson distribution and log link. For an estimation of our ability to detect actual effects, please see the power analysis for P3 below.
```{r expbol, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#UNREGISTERED ANALYSIS: first check whether first latency to land on the ground per video correlates between exploration environment and exploration object. The field season is now running in Woodland, CA, but the boldness data from Arizona is not yet done being video coded. The field team needs to know whether they can eliminate one of the exploration tests (they are currently conducting environment and object) because these assays are taking up too much time. The variable with the most data in it is Latency to first landing on the ground.
dexp <- read.csv("https://raw.githubusercontent.com/corinalogan/grackles/master/Files/Preregistrations/g_exploration_data_aviaryexplorationmotpersist.csv", header = T, sep = ",", stringsAsFactors = F)
#examine only Time 1 because not looking at intraindividual reliability and only the novel condition
dexp <- dexp[!dexp$Time==2 & !dexp$Condition=="Familiar",]
dexp$LatencyFirstLand <- as.numeric(dexp$LatencyFirstLand)
dexp$Test <- as.factor(dexp$Test)
dexp$BirdID <- as.factor(dexp$BirdID)
dexp$Time <- as.factor(dexp$Time)
#plot to see what the data show. Omit facet_wrap if you want to see all birds in one plot. Looks like there are differences between environment and object latencies for 12/19 birds. Generally faster to approach the environment than the object
library(ggplot2)
ggplot(dexp, aes(x = Test, y = LatencyFirstLand, colour = BirdID)) +
geom_point() +
facet_wrap( ~ BirdID)
library(lme4)
m.eo <- glmer(LatencyFirstLand ~ Test + (1|BirdID), family=gaussian, data=dexp[-which(is.na(dexp$LatencyFirstLand)),])
summary(m.eo) #negative correlation between environment and object latencies,r=-0.52
#However, if one of the exploration tests is repeatable and the other is not, this could explain the lack of correlation between the two conditions and tell us which test we want to continue with
### REGISTERED ANALYSES
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(LatencyBoldness ~ TimeOutsideNovelEnv + LatencyTableExpObject, family=poisson, data=persist2), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#GLM
eb <- glm(LatencyBoldness ~ TimeOutsideNovelEnv + LatencyTableExpObject, family="poisson", data=expbol)
#summary(eb)
eb2 <- summary(eb)
library(xtable)
eb2.table <- xtable(eb2)
library(knitr)
kable(eb2.table, caption="Table 2: Model selection output.", format="html", digits=2)
#Model Validation
library(MuMIn)
options(na.action = "na.fail")
base1 <- dredge(glm(NumberIncorrectTrialsReversal ~ AvgFunctionalTouchesExploration * Condition, family="poisson", data=persist2))
library(knitr)
kable(base1, caption="Table 3: Model selection output.")
```
**Model validation:** Determine whether the test model results are likely to be reliable given the data [@burnham2003model]. Compare Akaike weights (range: 0–1, the sum of all model weights equals 1; Akaike, 1981) between the test model and a base model (number of trials to reverse as the response variable and 1 as the explanatory variable) using the dredge function in the MuMIn package [@bates2012lme4]. If the best fitting model has a high Akaike weight (>0.89; [@burnham2003model]), then it indicates that the results are likely given the data. The Akaike weights indicate the best fitting model is the [base/test *- delete as appropriate*] model (Table 2).
#### *H1: P3: correlations between persistence measures*
**Analysis:** Generalized Linear Model (GLM; glm function, stats package) with a Poisson distribution and log link.
To determine our ability to detect actual effects, we ran a power analysis in G&ast;Power with the following settings: test family=F tests, statistical test=linear multiple regression: Fixed model (R^2 deviation from zero), type of power analysis=a priori, alpha error probability=0.05. We reduced the power to 0.70 and increased the effect size until the total sample size in the output matched our projected sample size (n=32). The protocol of the power analysis is here:
*Input:*
Effect size f² = 0,27
α err prob = 0,05
Power (1-β err prob - note: β=probability of making a Type II error) = 0,7
Number of predictors = 2
*Output:*
Noncentrality parameter λ = 8,6400000
Critical F = 3,3276545
Numerator df = 2
Denominator df = 29
Total sample size = 32
Actual power = 0,7047420
This means that, with our sample size of 32, we have a 70% chance of detecting a medium (approximated at f^2=0.15 by @cohen1988statistical) to large effect (approximated at f^2=0.35 by @cohen1988statistical).
```{r persist2, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
persist2 <- read.csv ("/Users/corina/GTGR/data/data_reverselatency.csv", header=T, sep=",", stringsAsFactors=F)
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(NumberIncorrectTrialsReversal ~ AvgFunctionalTouchesExploration * Condition, family=poisson, data=persist2), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#GLM
p2 <- glm(NumberIncorrectTrialsReversal ~ AvgFunctionalTouchesExploration * Condition, family="poisson", data=persist2)
#summary(p2)
sp2 <- summary(p2)
library(xtable)
sp2.table <- xtable(sp2)
library(knitr)
kable(sp2.table, caption="Table 2: Model selection output.", format="html", digits=2)
#Model Validation
library(MuMIn)
options(na.action = "na.fail")
base1 <- dredge(glm(NumberIncorrectTrialsReversal ~ AvgFunctionalTouchesExploration * Condition, family="poisson", data=persist2))
library(knitr)
kable(base1, caption="Table 3: Model selection output.")
```
**Model validation:** Determine whether the test model results are likely to be reliable given the data [@burnham2003model]. Compare Akaike weights (range: 0–1, the sum of all model weights equals 1; Akaike, 1981) between the test model and a base model (number of trials to reverse as the response variable and 1 as the explanatory variable) using the dredge function in the MuMIn package [@bates2012lme4]. If the best fitting model has a high Akaike weight (>0.89; [@burnham2003model]), then it indicates that the results are likely given the data. The Akaike weights indicate the best fitting model is the [base/test *- delete as appropriate*] model (Table 2).
#### *H2: P6: captive vs wild*
A GLMM (as in the repeatability analysis) will be conducted.
```{r capwildbold, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Boldness
cw <- read.csv ("/Users/corina/GTGR/data/data_individual_differences.csv", header=T, sep=",", stringsAsFactors=F)
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(LatencyBoldness ~ Context + AssayAttempts + (1|ID), family=poisson, data=cw), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#GLMM
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
bold <- MCMCglmm(LatencyBoldness ~ Context + AssayAttempts, random=~ID, family="poisson", data=cw, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(bold)
#autocorr(bold$Sol) #Did fixed effects converge?
#autocorr(bold$VCV) #Did random effects converge?
```
```{r capwildpers, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Boldness
cw <- read.csv ("/Users/corina/GTGR/data/data_individual_differences.csv", header=T, sep=",", stringsAsFactors=F)
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(NumberTouches ~ Context + AssayAttempts + (1|ID), family=poisson, data=cw), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#GLMM
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
bold <- MCMCglmm(NumberTouches ~ Context + AssayAttempts, random=~ID, family="poisson", data=cw, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(bold)
#autocorr(bold$Sol) #Did fixed effects converge?
#autocorr(bold$VCV) #Did random effects converge?
```
```{r capwildexp, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Exploration: novel environment
cw <- read.csv ("/Users/corina/GTGR/data/data_individual_differences.csv", header=T, sep=",", stringsAsFactors=F)
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(LatencyExpEnv ~ Context + AssayAttempts + (1|ID), family=poisson, data=cw), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#GLMM
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
bold <- MCMCglmm(LatencyExpEnv ~ Context + AssayAttempts, random=~ID, family="poisson", data=cw, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(bold)
#autocorr(bold$Sol) #Did fixed effects converge?
#autocorr(bold$VCV) #Did random effects converge?
#Exploration: novel object
cw <- read.csv ("/Users/corina/GTGR/data/data_individual_differences.csv", header=T, sep=",", stringsAsFactors=F)
# DATA CHECKING
library(DHARMa)
library(lme4)
simulationOutput <- simulateResiduals(fittedModel = glmer(LatencyTableExpObject ~ Context + AssayAttempts + (1|ID), family=poisson, data=cw), n=250) #250 simulations, but if want higher precision change n>1000
simulationOutput$scaledResiduals #Expect a flat distribution of the overall residuals, and uniformity in y direction if plotted against any predictor
testDispersion(simulationOutput) #if under- or over-dispersed, then p-value<0.05, but then check the dispersion parameter and try to determine what in the model could be the cause and address it there, also check for zero inflation
testZeroInflation(simulationOutput) #compare expected vs observed zeros, not zero-inflated if p<0.05
testUniformity(simulationOutput) #check for heteroscedasticity ("a systematic dependency of the dispersion / variance on another variable in the model" Hartig, https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html), which is indicated if dots aren't on the red line and p<0.05. Also...
plot(simulationOutput) #...there should be no pattern in the data points in the right panel
plotResiduals(ReverseNumber, simulationOutput$scaledResiduals) #plot the residuals against other predictors (in cases when there is more than 1 fixed effect) - can't get this code to work yet
#GLMM
library(MCMCglmm)
prior = list(R=list(R1=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0)))
bold <- MCMCglmm(LatencyTableExpObject ~ Context + AssayAttempts, random=~ID, family="poisson", data=cw, verbose=F, prior=prior, nitt=13000, thin=10, burnin=3000)
summary(bold)
#autocorr(bold$Sol) #Did fixed effects converge?
#autocorr(bold$VCV) #Did random effects converge?
```
#### *Alternative Analyses*
We anticipate that we will want to run additional/different analyses after reading @statrethinkingbook. We will revise this preregistration to include these new analyses before conducting the analyses above.
### F. ETHICS
This research is carried out in accordance with permits from the:
1) US Fish and Wildlife Service (scientific collecting permit number MB76700A-0,1,2)
2) US Geological Survey Bird Banding Laboratory (federal bird banding permit number 23872)
3) Arizona Game and Fish Department (scientific collecting license number SP594338 [2017] and SP606267 [2018])
4) Institutional Animal Care and Use Committee at Arizona State University (protocol number 17-1594R)
5) University of Cambridge ethical review process (non-regulated use of animals in scientific procedures: zoo4/17)
### G. AUTHOR CONTRIBUTIONS
**McCune:** Hypothesis development, data collection, data analysis and interpretation, write up, revising/editing.
**MacPherson:** Data collection, data interpretation, revising/editing.
**Rowney:** Data collection, data interpretation, revising/editing.
**Bergeron:** Data collection, data interpretation, revising/editing.
**Folsom:** Data collection, data interpretation, revising/editing.
**Deffner:** Data analysis (Flexibility comprehensive model), data interpretation, revising/editing.
**Breen:** Data collection (video coding), data interpretation, revising/editing.
**Logan:** Hypothesis development, data collection, data analysis and interpretation, revising/editing, materials/funding.
### H. FUNDING
This research is funded by the Department of Human Behavior, Ecology and Culture at the Max Planck Institute for Evolutionary Anthropology, and by a Leverhulme Early Career Research Fellowship to Logan in 2017-2018.
###I. CONFLICT OF INTEREST DISCLOSURE
We, the authors, declare that we have no financial conflicts of interest with the content of this article. Corina Logan is a Recommender and on the Managing Board at PCI Ecology.
### J. ACKNOWLEDGEMENTS
We thank Dieter Lukas for help polishing the predictions; Ben Trumble for providing us with a wet lab at Arizona State University and Angela Bond for lab support; Melissa Wilson for sponsoring our affiliations at Arizona State University and lending lab equipment; Kevin Langergraber for serving as local PI on the ASU IACUC; Kristine Johnson for technical advice on great-tailed grackles; Arizona State University School of Life Sciences Department Animal Care and Technologies for providing space for our aviaries and for their excellent support of our daily activities; Julia Cissewski for tirelessly solving problems involving financial transactions and contracts; Richard McElreath for project support; Aaron Blackwell and Ken Kosik for being the UCSB sponsors of the Cooperation Agreement with the Max Planck Institute for Evolutionary Anthropology; Jeremy Van Cleve, our Recommender at PCI Ecology, and two anonymous reviewers for their wonderful feedback; Vincent Kiepsch and Sierra Planck for interobserver reliability video coding: Sawyer Lung for field support; and our research assistants: Aelin Mayer, Nancy Rodriguez, Brianna Thomas, Aldora Messinger, Elysia Mamola, Michael Guillen, Rita Barakat, Adriana Boderash, Olateju Ojekunle, August Sevchik, Justin Huynh, Jennifer Berens, Amanda Overholt, Michael Pickett, Sam Munoz, Sam Bowser, Emily Blackwell, Kaylee Delcid, Sofija Savic, Brynna Hood, Sierra Planck, and Elise Lange.
### K. [REFERENCES](MyLibrary.bib)