Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
770 lines (546 sloc) 75.5 KB
---
title: Is behavioral flexibility related to foraging and social behavior in a rapidly expanding species?
author: '[Dr. Corina Logan](http://CorinaLogan.com) (Max Planck Institute for Evolutionary Anthropology, MPI EVA, corina_logan@eva.mpg.de), Luisa Bergeron (University of California Santa Barbara, UCSB), Dr. Kelsey McCune (UCSB), Melissa Folsom (MPI EVA), [Dr. Dieter Lukas](http://dieterlukas.strikingly.com) (MPI EVA)'
date: '`r Sys.Date()`'
output:
md_document:
toc: true
html_document:
toc: true
toc_depth: 4
toc_float:
collapsed: false
code_folding: hide
pdf_document:
keep_tex: yes
latex_engine: xelatex
github_document:
toc: true
bibliography: /Users/corina/GitHub/grackles/Files/MyLibrary.bib
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r}
#Make code wrap text so it doesn't go off the page when Knitting to PDF
library(knitr)
opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=TRUE)
```
***Click [here](https://github.com/corinalogan/grackles/blob/master/README.md) to navigate to the version-tracked reproducible manuscript (.Rmd file)***
###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 (e.g., behavior, invasion success, diet generalism, foraging techniques, foraging innovations, mortality, brain size), 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 (i.e., they increased their nesting range by over 5500% between 1880 and 2000 [@wehtje2003range], [@peer2011invasion]) (Fig. 1). Foraging behavior is considered central to the rapid geographic range expansion of this species and it is thought that they have been so successful by following human urban and agricultural corridors ([@wehtje2003range], [@peer2011invasion]). Therefore, as humans continue to modify landscapes, this increases the amount of suitable grackle habitat. We expect this species to be behaviorally flexible because they are fast at reversal learning (@logan2016behavioral), they often encounter human-made "puzzle boxes" in the wild as they attempt to open packaging to access food when digging through garbage cans and eating at outdoor cafes, and they may track resources across time and space. Results will allow us to determine whether, as predicted by hypotheses and cross-species correlational data, in this expanding species, individual-level variation in flexibility is linked with diet breadth, foraging proficiency, social interactions, habitat use, and movement into new geographic areas. **This investigation**: In this piece of the long-term project, we will assess whether individual performance in experiments that assess behavioral flexibility relates to individual variation in ecological and social behavior in the natural environment. In particular, we aim to determine whether the more behaviorally flexible (measured by reversal learning and solution switching on a multi-access box in a separate [preregistration](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md)) grackles have more flexible foraging behavior (eat a larger number of different foods, use a wider variety of foraging techniques), are more flexible in their habitat use (are found in more diverse habitat types, disperse farther from their natal area), and are more flexible in their social relationships (have more or stronger social bonds particularly with less related individuals). We will be able to compare the grackle's ability to adapt behavior according to social context with data from other species, as well as determine whether it is linked with measures of flexibility in asocial contexts.
![Figure 1. Five-year project overview.](https://github.com/corinalogan/grackles/blob/master/Files/GrackleProjectTimeline.png)
###A. STATE OF THE DATA
This preregistration was written (Dec 2017) prior to collecting any data, and submitted to PCI Ecology for peer review (Oct 2018) after 46 practice focal follows had been conducted (we will not use this data in our analyses below because they were for learning how to do focals and for testing interobserver reliability). We received the peer reviews (Mar 2019, after 3 more practice follows had been conducted (these data will not be used in analyses) and 13 real follows had been conducted), revised, and resubmitted to PCI Ecology (Apr 2019, 3 more real follows had been conducted) after 16 focal follows had been conducted.
###B. PARTITIONING THE RESULTS
We may decide to present the results from different hypotheses in separate papers.
###C. HYPOTHESIS
#### H1: [Behavioral flexibility](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md) (see @mikhalevich_is_2017 for a detailed definition; measured by reversal learning [where they must learn to prefer one of two options that contain food and then reverse this preference] and switching between options on a multi-access box (where they must learn to switch to a new option, out of four available options, when an option becomes non-functional) in aviaries) is related to foraging behavior (measured with focal follows using this [ethogram](https://docs.google.com/spreadsheets/d/1N8wsA3geaRGlMjRxYTRpdG2i5oCXNGq9zBlTnj02Gho/edit?usp=sharing)) in wild individuals (after their release from the aviaries). We expect this species to be behaviorally flexible because they are fast at reversal learning [@logan2016behavioral], they often encounter human-made "puzzle boxes" in the wild as they attempt to open packaging to access food when digging through garbage cans and eating at outdoor cafes, and they may track resources across time and space. Foraging behavior is considered central to the rapid geographic range expansion of this species and it is thought that they have been so successful by following human urban and agricultural corridors ([@wehtje2003range], [@peer2011invasion]). Therefore, as humans continue to modify landscapes, this increases the amount of suitable grackle habitat.
**Prediction 1:** Individuals that are faster to reverse preferences on a reversal learning task and who also have lower latencies to switch to solving new loci after previously solved loci become unavailable (multi-access box) will eat a larger number of different foods and use a wider variety of foraging techniques in the wild, validating the cross-species correlational finding that technique breadth (@overington2009technical) and diet breadth (@ducatez2015ecological) is associated with flexibility.
**P1 alternative 1:** If there is no correlation, this suggests that flexibility as we measured it represents a trait that is not related to the number of foods eaten and foraging techniques used. Flexibility may not necessarily be associated with diet and foraging technique breadth because flexibility could be constrained in a foraging context due to social competition (e.g., subordinates are outcompeted while foraging and thus try new foods and techniques) or ecological limitations (e.g., constrained by what is available). Additional research would be required to determine the factors that might constrain foraging behavior.
**P1 alternative 2:** If there is a negative correlation between flexibility and the number of different foods eaten, this might indicate that the more flexible individuals target particular food items.
**P1 alternative 3:** If there is a negative correlation between flexibility and the number of foraging techniques, this could indicate that the more flexible individuals use particular (potentially more effective) techniques.
**P2:** Individuals whose [flexibility](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md) has been increased experimentally will consume a larger number of foods and use more foraging techniques (measured with focal follows) than individuals whose flexibility has not been manipulated. This would further validate that flexibility is related to diet breadth and foraging techniques.
**P2 alternative 1:** If the flexibility manipulation does not work in that those individuals in the experimental condition do not decrease their reversal learning speeds more than control individuals, then we will rely on the general individual variation in flexibility and how it relates to foraging in the wild (as in P1).
**P3:** The proportion of their diet that is human foods and the proportion of their foraging techniques involving human foods is higher for the more flexible individuals, who will consistently occur in locations closer to known outdoor human food locations like picnic areas and outdoor cafe seating (measured as the repeatability of the individual's distance from cafes across multiple separate focal follows). For the diet, this is potentially due to A) having stayed in their parent's home range (i.e., they eat human food because it happens to be more prevalent in their home range than in other home ranges; local specialization) or B) because these individuals move around to seek out such opportunities (potentially seeking out habitat edges within their population). For the foraging techniques, this is potentially due to human foods and their packaging changing at a faster rate than natural foods and prey items and their accessibility. Foods eaten and foraging techniques used will be recorded during focal follows. Because this species is highly associated with human-modified landscapes, it is likely that consuming human foods is part of the reason for this association, and that flexible individuals are better at solving these human-made "puzzle boxes" to access food.
**P3 alternative 1:** There is no correlation between an individual's flexibility and the proportion of human foods in their diet, potentially because A) their daily range sizes encompass many different food resources, including human foods (though they are likely not specialized on human foods), and B) some less flexible individuals might specialize on human foods.
**P3 alternative 2:** There is a negative correlation between an individual's flexibility and the proportion of human foods in their diet, potentially because some of the less flexible individuals might specialize on human foods, thus increasing their consumption above that of the more flexible individuals.
#### H2: [Behavioral flexibility](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md) (see @mikhalevich_is_2017 for a detailed definition; measured by reversal learning and switching between options on a multi-access box in aviaries) is related to social behavior (measured with focal follows using this [ethogram](https://docs.google.com/spreadsheets/d/1N8wsA3geaRGlMjRxYTRpdG2i5oCXNGq9zBlTnj02Gho/edit?usp=sharing)) in wild individuals. To give an example of the types of social relationships this species engages in, in terms of male social relationships, @johnson2000male found during the breeding season in a population in Texas that one or more territorial males defend a territory with several nests from females, that non-territory holding resident males will queue to gain access to a territory, and that transient males move from colony to colony. There could be varying needs for males to manage their relationships with each other and flexibility could potentially play a role in such management.
**P4** Flexible individuals are more likely to have a greater number of bonds OR stronger bonds with others, in particular with individuals who are less related, potentially because they are better able to adjust their behavior to that of an affiliate. Social bonds are measured using the focal follow method to sample affiliative and aggressive behaviors.
**P4 alternative 1:** Individual flexibility is not related to the number or strength of social bonds, potentially because all individuals are able to form bonds with like individuals, including the less flexible individuals.
**P4 alternative 2:** Flexible individuals may have fewer affiliates or be less likely to regularly affiliate with the same individuals, potentially because they frequently change their behavior and are difficult to associate with. We are not able to test this alternative in this study, but could propose experimental designs for future research if this alternative is supported by the data.
#### H3: [Behavioral flexibility](https://github.com/corinalogan/grackles/blob/master/EasyToReadFiles/g_flexmanip.md) (see @mikhalevich_is_2017 for a detailed definition; measured by reversal learning and switching between options on a multi-access box in aviaries) will differ in their use of microhabitats within human-modified landscapes (substrate qualification during each focal follow), but the macrohabitat (square kilometer) of each population will not differ in human population density (measured with a GPS point for each focal follow after their release from the aviaries; we measure microhabitat types according to the primary substrate the focal follow occurred on: grass, gravel, tree, cement cafe, dumpster, [additional substrates will be added as they are encountered]). Although we were only able to find this species in association with human-modified landscapes based on eBird sightings (i.e., there appear to be no forest-based populations), individuals could use these landscapes in a variety of ways. For example they could specialize on particular foods or at particular types of locations (e.g., foraging exclusively at cafes or in grassy areas), they could generalize across all foods and location types, or they might fall somewhere in between these extremes.
**P5:** Individuals immigrating into a population are more likely to be flexible, potentially because they need to learn how to obtain resources in an unfamiliar area. Immigrants are individuals who carry many genetic variants (identified using ddRADseq) that are not found in other individuals in this population.
**P5 alternative:** Individuals immigrating into a population are not more likely to be flexible, potentially because the human urban environment is comparable across landscapes.
**P6:** Flexible individuals will be found more regularly in a wider diversity microhabitats (human-modified substrates including cement, dumpsters or cafes; or natural substrates including grass, shrubs, and trees [additional substrates will be added as they are encountered]) during focal follows.
**P6 alternative:** Flexibility is not associated with presence in diverse microhabitats because the more flexible individuals might specialize in specific foraging strategies.
**P7:** There will be no difference in human population density among the sites for our three grackle populations because all great-tailed grackle populations are highly associated with human-modified landscapes. Human population density per square mile data will be obtained from census information (US census bureau: https://www.census.gov/quickfacts/fact/note/US/LND110210, still looking for a source for Central American countries)
###D. METHODS
####**Planned sample**
Great-tailed grackles (n > 200) will be caught in the wild at three field sites across their geographic range: the center of their original range (at a site to be determined in Central America), the middle of the northward expanding edge (Tempe, Arizona USA), and on the northern expanding edge (to be determined). Individuals will be identified using colored leg bands in unique combinations, their data collected (blood, feathers, and biometrics), and then they will be released back to the wild. Some individuals (64-100) will be brought temporarily into aviaries for behavioral testing, and then they will be released back to the wild where the data for this study will be collected.
**Sample size rationale**
We will band as many birds as possible throughout the year, and conduct behavioral tests in aviaries (during the non-breeding seasons; approximately September through March) and in the wild (year-round) on as many birds as we can at each field site. The minimum sample size will be 200 banded birds and 60 behavior-tested birds in total, however we expect to be able to sample many more.
**Data collection stopping rule**
We will stop collecting data in April 2023 when the current funding ends, or after data have been collected at all three field sites, whichever date comes first.
####**Open materials**
[Ethogram](https://docs.google.com/spreadsheets/d/1N8wsA3geaRGlMjRxYTRpdG2i5oCXNGq9zBlTnj02Gho/edit?usp=sharing) for Prim8.
[Individuals](https://docs.google.com/spreadsheets/d/1Lr0pwsmdnpVM8X2Fyoj9EIGa3zOY1WCZlntW7e0Ui_Y/edit?usp=sharing) for Prim8.
####**Open data**
When the study is complete, the data will be published in the Knowledge Network for Biocomplexity's data repository.
####**Randomization and counterbalancing**
No randomization or counterbalancing is involved in this study.
####**Blinding of conditions during analysis**
No blinding is involved in this study.
####**Dependent variables**
***P1-P2***
1) Number of different foods eaten in the first X minutes (X=the sum of the total observation time per individual, using the individual who had the lowest sum to equalize observation time across individuals)
2) Number of foraging techniques used (based on Table 1 in @overington2009technical) in the first X minutes (X=the sum of the total observation time per individual, using the individual who had the lowest sum to equalize observation time across individuals)
One model will be run per dependent variable.
***P3: flexible = more human foods***
1) Proportion of diet per individual that is human food
***P4: flexible = a greater number of bonds or stronger bonds***
1) Strength of the maximum bond (calculated as the half-weight index based on association behavior during focal follows).
2) Individual strength (the sum of all bonds an individual has; @wey2008social)
3) Individual degree (maximum number of other individuals that the focal subject associated with; @wey2008social).
4) Male shares territory with another male: yes, no
5) Relatedness for the strongest bond (measured following the protocol in @thrasher2018double to estimate pairwise relatedness between all individuals based on the extent of sharing of genetic variants as determined by ddRADseq)
***P5: flexible = immigrants***
1) Probability of being an immigrant (measured using maximum-likelihood-based individual assignment probabilities calculated using an admixture program with the ddRADseq data)
***P6: flexible = wider range of habitats***
1) Evenness in the proportion of time spent in each habitat type (grass, gravel and other natural substrate, cement, cafe, dumpster)
####**Independent variables**
***P1-P4 and P6***
1) Flexibility 1: **Number of trials to reverse** a preference in the last reversal (in the reversal learning experiment) an individual experienced (individuals in the flexibility control group only experience 1 reversal so this data will come from their first and only reversal; individuals in the flexibility manipulation group experience serial reversals until they pass a certain criterion, therefore we will only use data from their most recent reversal). 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).
2) Flexibility 3: If the number of trials to reverse a preference does not positively correlate with the number of trials to attempt or solve (meet criterion) new loci on the multi-access box (an additional measure of behavioral flexibility), then the **number of trials to solve** and the **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 4: 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. It will be based on a Bayesian estimate of the reduction in error across trials estimated from the number of correct choices from the beginning of each reversal. 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 flexibility measures 1 through 3. If this ends up being the case, we will modify the code in the analysis plan below to reflect this change.
4) Flexibility manipulation: control, manipulated
5) Dominance rank
6) Condition: aviary-tested, not-aviary-tested
7) ID (random effect because multiple measures per individual)
8) Population: center (Central America), middle (Arizona), edge (northern US) (random effect because each population might have a different slope)
***P5***
1) Flexibility 1: **Number of trials to reverse** a preference in the last reversal an individual experienced (reversal learning; 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).
2) Flexibility 3: If the number of trials to reverse a preference does not positively correlate with the latency to attempt or solve new loci on the multi-access box (an additional measure of behavioral flexibility), then the **average latency to solve** and the **average latency 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 4: A Bayesian 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 flexibility measures 1 through 3. If this ends up being the case, we will modify the code in the analysis plan below to reflect this change.
4) Population: center (Central America), middle (Arizona), edge (northern US) (random effect because each population might have a different slope)
####**Focal follow interobserver reliability**
To be able to conduct focal follows (methods as in @altmann1974observational), a coder must pass interobserver reliability before the data they collect is used in the data set. To pass, coders must have an intra-class correlation (ICC; @hutcheon2010random) of 0.90 or greater based on at least six 10-min focal follows where both coders recorded the behavior of the same focal individual at the same time.
Luisa Bergeron was the first person to conduct focal follows, therefore she trained Kelsey McCune and Melissa Folsom until they passed interobserver reliability (on 10 June 2019) for each of the 6 variables listed in the preregistrations.
```{r ior, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
ff <- read.csv ("/Users/corina/GTGR/data/data_flexforaging.csv", header=T, sep=",", stringsAsFactors=F)
#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."
#Download the datasheet as a .csv file and put in a folder where you know the file path
#Load the datasheet here
data <- read.csv ("/Users/corina/ownCloud/Documents/Grackles/Data/ASU/InterObsRelKel.csv", header = TRUE, sep = ",", stringsAsFactors = FALSE)
data #Check to make sure it looks right
#LOAD the irr package, and run it in R
library(irr)
### ICCs for agreement between the 2 coders for 6 variables
#Note: c(4,5) is telling R to look at columns 4 ("1NumberForagingTechniques") and 5 ("2NumberForagingTechniques") and compare them
#NumberDifferentFoodsEaten
icc(data[,c(2,3)], model="oneway", type="agreement", unit="single", conf.level=0.95)
#NumberDifferentForagingTechniques
icc(data[,c(4,5)], model="oneway", type="agreement", unit="single", conf.level=0.95)
#NumberAffiliativeInteractions
icc(data[,c(6,7)], model="oneway", type="agreement", unit="single", conf.level=0.95)
#NumberAggressiveInteractions
icc(data[,c(8,9)], model="oneway", type="agreement", unit="single", conf.level=0.95)
#NumberTimesSubjectInitiatedAggression
icc(data[,c(10,11)], model="oneway", type="agreement", unit="single", conf.level=0.95)
#Microhabitat
icc(data[,c(12,13)], model="oneway", type="agreement", unit="single", conf.level=0.95)
```
**ICCs for Melissa (n=6 focal follows, [data](https://github.com/corinalogan/grackles/blob/master/Files/Data/data_InterObsRelMel.csv)):**
Different Foods Eaten = 1.00
Diffent Foraging Techniques = 1.00
Number of Affiliative Interactions = 1.00
Number of Aggressive Interactions = 0.96
Number of Initiated Aggressive Interactions = 0.94
Microhabitat = 1.00
**ICCs for Kelsey (n=6 focal follows, [data](https://github.com/corinalogan/grackles/blob/master/Files/Data/data_InterObsRelKel.csv)):**
Different Foods Eaten = 1.00
Different Foraging Techniques = 0.97
Number of Affiliative Interactions = 0.96
Number of Aggressive Interactions = 1.00
Number of Initiated Aggressive Interactions = 1.00
Microhabitat = 1.00
NOTE: the ICCs for the variable Different Foods Eaten for these focal follows was originally 0.63 (Melissa) and 0.64 (Kelsey) because Melissa and Kelsey recorded a "bug" being eaten while Luisa recorded no food type because she couldn't identify it to a more specific category. At this point, we decided that we would prefer to enter a general category for food type rather than having no information about what was eaten. Therefore, this data point was removed from the interobserver reliability analysis. This resulted in ICCs of 1.00 for both Kelsey and Melissa on the Different Foods Eaten variable because they matched Luisa in the other food type data points.
###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.
We will analyze data for females and males separately because each sex has a distinct natural history.
####*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 will revise the Analysis Plan after reading @statrethinkingbook. We currently don’t have any estimates for any of our measured variables because these tests have never been done in grackles and we have not encountered previous research that has manipulated flexibility in this way. We will use Bayesian analyses to estimate our likely confidence in the results given simulated data. We will revise this preregistration to include these new analyses before conducting the planned analyses on our actual data. Based on the simulations, we might adapt the number of focal follows per individual or decide to collect much more data just with the aviary-tested birds to increase the amount of information per individual.
####*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].
####*P1-P2*
**Analysis:** 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; [@hadfield2010mcmc]) with a Poisson distribution and log link using 130,000 iterations with a thinning interval of 10, a burnin of 30,000, and minimal priors (V=1, nu=0) [@hadfield2014coursenotes]. We will ensure the GLMM shows acceptable convergence (lag time autocorrelation values <0.01; [@hadfield2010mcmc]), and adjust parameters if necessary to meet this criterion. We will determine whether an independent variable had an effect or not using the Estimate in the full model.
```{r p1p2, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
ff <- read.csv ("/Users/corina/GTGR/data/data_flexforaging.csv", header=T, sep=",", stringsAsFactors=F)
#Separate the sexes
fem <- ff[ff$Sex=="f",]
mal <- ff[ff$Sex=="m",]
#Factor the random effect variables
ID <- as.factor(fem$ID)
Population <- as.factor(fem$Population)
ID <- as.factor(mal$ID)
Population <- as.factor(mal$Population)
# DATA CHECKING
library(DHARMa)
library(lme4)
#Data checking for GLMM 1 females
simulationOutput <- simulateResiduals(fittedModel = glmer(NumberFoodsEaten ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + (1|ID) + (1|Population), family=poisson, data=fem), 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(FlexRatio, 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
#Data checking for GLMM 1 males
simulationOutput <- simulateResiduals(fittedModel = glmer(NumberFoodsEaten ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + (1|ID) + (1|Population), family=poisson, data=mal), 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(FlexRatio, 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
#Data checking for GLMM 2 females
simulationOutput <- simulateResiduals(fittedModel = glmer(NumberForagingTechniques ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=fem), 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(FlexRatio, 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
#Data checking for GLMM 2 males
simulationOutput <- simulateResiduals(fittedModel = glmer(NumberForagingTechniques ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=mal), 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(FlexRatio, 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),R2=list(V=1,nu=0),R3=list(V=1,nu=0),R4=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0),G2=list(V=1,nu=0)))
#GLMM 1 with response variable = NumberFoodsEaten
#females
f1 <- MCMCglmm(NumberFoodsEaten ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f1)
#autocorr(f1$Sol) #Did fixed effects converge?
#autocorr(f1$VCV) #Did random effects converge?
#males
f2 <- MCMCglmm(NumberFoodsEaten ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f2)
#autocorr(f2$Sol) #Did fixed effects converge?
#autocorr(f2$VCV) #Did random effects converge?
#GLMM 2 with response variable = NumberForagingTechniques
#female
f3 <- MCMCglmm(NumberForagingTechniques ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f3)
#autocorr(f3$Sol) #Did fixed effects converge?
#autocorr(f3$VCV) #Did random effects converge?
#male
f4 <- MCMCglmm(NumberForagingTechniques ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f4)
#autocorr(f4$Sol) #Did fixed effects converge?
#autocorr(f4$VCV) #Did random effects converge?
```
We will quantify the number of different food types and foraging techniques during focal follows according to the [ethogram](https://docs.google.com/spreadsheets/d/1N8wsA3geaRGlMjRxYTRpdG2i5oCXNGq9zBlTnj02Gho/edit?usp=sharing). If a grackle forages during a focal follow we record WHAT it eats, as well as HOW the bird is searching for food. Foraging techniques include: Flipping over objects (flip), digging in ground with bill or feet (dig), sweeping head back and forth [i.e., actually sweeping the bill across the substrate; (sw)], extracting from a substrate (ex), lowers body posture to be parallel to ground to stalk/catch prey from air, from ground, from tree, etc (sc), using gaping bill to search through substrate (ga), lifting or nudging objects with bill (ln).
####*P3: flexible = more human foods*
**Analysis:** A GLMM was conducted as in P1-P2, except this GLMM used a binomial distribution (called "categorical" in MCMCglmm) due to the response variable being a proportion.
```{r p3, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
ff <- read.csv ("/Users/corina/GTGR/data/data_flexforaging.csv", header=T, sep=",", stringsAsFactors=F)
#Separate the sexes
fem <- ff[ff$Sex=="f",]
mal <- ff[ff$Sex=="m",]
#Factor the random effect variables
ID <- as.factor(fem$ID)
Population <- as.factor(fem$Population)
ID <- as.factor(mal$ID)
Population <- as.factor(mal$Population)
# DATA CHECKING
library(DHARMa)
library(lme4)
#Data checking for female GLMM
simulationOutput <- simulateResiduals(fittedModel = glmer(ProportionFoodsEatenHumanFood ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=binomial, data=fem), 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(ProportionFoodsEatenHumanFood, 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
#Data checking for male GLMM
simulationOutput <- simulateResiduals(fittedModel = glmer(ProportionFoodsEatenHumanFood ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=binomial, data=mal), 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(ProportionFoodsEatenHumanFood, 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),R2=list(V=1,nu=0),R3=list(V=1,nu=0),R4=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0),G2=list(V=1,nu=0)))
#GLMM with response variable = NumberFoodsEaten
#females
f1 <- MCMCglmm(ProportionFoodsEatenHumanFood ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="categorical", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f1)
#autocorr(f1$Sol) #Did fixed effects converge?
#autocorr(f1$VCV) #Did random effects converge?
#males
f2 <- MCMCglmm(ProportionFoodsEatenHumanFood ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="categorical", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f2)
#autocorr(f2$Sol) #Did fixed effects converge?
#autocorr(f2$VCV) #Did random effects converge?
# Is distance from outdoor cafes a repeatable trait within individuals, as measured with a GPS point of individual locations during several separate focal follows
d1 <- rpt(dist ~ Sex + (1 | ID), grname = "ID", data = ff, datatype = "poisson", nboot = 1000, npermut = 300)
summary(d1)
```
####*P4: flexible = a greater number of bonds or stronger bonds*
**Analysis:** A GLMM was conducted as in P1-P2.
To quantify social relationships, we will conduct at least four 10-minute focal follows on each subject. We will find subjects in the wild because we will attach radio transmitter tags to all grackles that are released from the aviaries upon completion of their test battery. To ensure we fully sample social and foraging behavior, we will prioritize conducting focal follows on these tagged grackles for which we have a much larger amount of individualized data, including multiple measures of flexibility. We will also sample many other banded grackles that were never tested in the aviaries, and thus do not have measures of flexibility. By conducting focal follows on grackles that were not in captivity, we can verify that the time in captivity had no effect on grackle social behavior after release because aviary-tested birds should be indistinguishable from non-aviary-tested birds in these analyses.
To measure affiliative bonds, during each focal follow we will record when another grackle comes within one body length of the focal bird (and does not engage in aggressive interactions). In case we do not observe enough of these close associations, we will also record when another grackle comes within 3m of the focal subject (and does not engage in aggressive interactions). Finally, we will conduct a scan sample at the end of the follow to determine group size as the number of other grackles within 10m of the focal individual. Unbanded grackles that are seen in proximity of the focal individual will be recorded and included in the count of group size and individual degree (the number of unique associates), but because we cannot distinguish unbanded individuals from each other, we will exclude unbanded bird data from calculations of an individual’s summed bond strengths (see details in the next paragraph).
We will also measure aggressive behavioral interactions, as indicated in our ethogram. The outcome of these dyadic interactions will be used to create our index of dominance (wins - losses / wins + losses).
We will conduct subsequent follows on the same individual only when 3 or more weeks have passed since the previous focal follow to prevent temporal autocorrelation in behavior (@whitehead2008analyzing). From the data sheet of dyadic associations during focal follows, we will create a matrix of association strengths between all banded grackles by calculating the Half-Weight association index. This index determines association strength based on the proportion of observations in which two individuals are seen together versus separately, and accounts for bias arising from subjects that are more likely to be observed separately rather than together in the same group (@cairns1987comparison). From the matrix of association values, we will use the R package igraph (@csardi2006igraph) to create a social network, and calculate each individual’s strength (sum of all association values) and degree (maximum number of unique associates) values (@croft2008exploring).
Social network data are not independent (@croft2011hypothesis), therefore, to determine whether individuals are associating non-randomly based on flexibility (i.e., association strength between two grackles is larger than would be expected by random chance), we will compare our model results to those obtained from random networks. To make a random network, we will use the R package sna (@butts2016sna) and the function “rperm” to randomly rearrange the association strengths (edges) of grackles in our network. We will conduct this edge randomization (called permutation) 10,000 times to create our sample of random networks. We will then re-calculate our dependent variables from the random networks and re-run the same models (as in @croft2011hypothesis and @whitehead2005testing). We will conclude that social bonds are significantly related to flexibility if the coefficients describing the relationship in our observed data are in the top 2.5% and bottom 2.5% of the coefficients resulting from models run on the random networks.
```{r p4, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
ff <- read.csv ("/Users/corina/GTGR/data/data_flexforaging.csv", header=T, sep=",", stringsAsFactors=F)
#Separate the sexes
fem <- ff[ff$Sex=="f",]
mal <- ff[ff$Sex=="m",]
#Factor the random effect variables
ID <- as.factor(fem$ID)
Population <- as.factor(fem$Population)
ID <- as.factor(mal$ID)
Population <- as.factor(mal$Population)
# DATA CHECKING
library(DHARMa)
library(lme4)
#Data checking for female GLMM 1
simulationOutput <- simulateResiduals(fittedModel = glmer(MaxBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=fem), 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(MaxBondStrength, 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
#Data checking for male GLMM 1
simulationOutput <- simulateResiduals(fittedModel = glmer(MaxBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=mal), 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(MaxBondStrength, 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
#Data checking for female GLMM 2
simulationOutput <- simulateResiduals(fittedModel = glmer(SumBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=fem), 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(SumBondStrength, 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
#Data checking for male GLMM 2
simulationOutput <- simulateResiduals(fittedModel = glmer(SumBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=mal), 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(SumBondStrength, 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
#Data checking for female GLMM 3
simulationOutput <- simulateResiduals(fittedModel = glmer(Degree ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=fem), 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(Degree, 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
#Data checking for male GLMM 3
simulationOutput <- simulateResiduals(fittedModel = glmer(Degree ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=mal), 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(Degree, 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
#Data checking for female GLMM 4
simulationOutput <- simulateResiduals(fittedModel = glmer(MaleSharesTerritory ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=fem), 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(MaleSharesTerritory, 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
#Data checking for male GLMM 4
simulationOutput <- simulateResiduals(fittedModel = glmer(MaleSharesTerritory ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=mal), 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(MaleSharesTerritory, 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
#Data checking for female GLMM 5
simulationOutput <- simulateResiduals(fittedModel = glmer(RelatednessForMaxBond ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=fem), 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(RelatednessForMaxBond, 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
#Data checking for male GLMM 5
simulationOutput <- simulateResiduals(fittedModel = glmer(RelatednessForMaxBond ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|ID) + (1|Population), family=poisson, data=mal), 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(RelatednessForMaxBond, 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),R2=list(V=1,nu=0),R3=list(V=1,nu=0),R4=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0),G2=list(V=1,nu=0)))
#GLMM 1 with response variable = MaxBondStrength
#females
f1 <- MCMCglmm(MaxBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f1)
#autocorr(f1$Sol) #Did fixed effects converge?
#autocorr(f1$VCV) #Did random effects converge?
#males
f2 <- MCMCglmm(MaxBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f2)
#autocorr(f2$Sol) #Did fixed effects converge?
#autocorr(f2$VCV) #Did random effects converge?
#GLMM 2 with response variable = AvgBondStrength
#female
f3 <- MCMCglmm(SumBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f3)
#autocorr(f3$Sol) #Did fixed effects converge?
#autocorr(f3$VCV) #Did random effects converge?
#male
f4 <- MCMCglmm(SumBondStrength ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f4)
#autocorr(f4$Sol) #Did fixed effects converge?
#autocorr(f4$VCV) #Did random effects converge?
#GLMM 3 with response variable = Degree
#female
f5 <- MCMCglmm(Degree ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f5)
#autocorr(f5$Sol) #Did fixed effects converge?
#autocorr(f5$VCV) #Did random effects converge?
#male
f6 <- MCMCglmm(Degree ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f6)
#autocorr(f6$Sol) #Did fixed effects converge?
#autocorr(f6$VCV) #Did random effects converge?
#GLMM 4 with response variable = MaleSharesTerritory
#female
f7 <- MCMCglmm(MaleSharesTerritory ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f7)
#autocorr(f7$Sol) #Did fixed effects converge?
#autocorr(f7$VCV) #Did random effects converge?
#male
f8 <- MCMCglmm(MaleSharesTerritory ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f8)
#autocorr(f8$Sol) #Did fixed effects converge?
#autocorr(f8$VCV) #Did random effects converge?
#GLMM 5 with response variable = RelatednessForMaxBond
#female
f9 <- MCMCglmm(RelatednessForMaxBond ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f9)
#autocorr(f9$Sol) #Did fixed effects converge?
#autocorr(f9$VCV) #Did random effects converge?
#male
f10 <- MCMCglmm(RelatednessForMaxBond ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f4)
#autocorr(f10$Sol) #Did fixed effects converge?
#autocorr(f10$VCV) #Did random effects converge?
```
####*P5: flexible = immigrants*
**Analysis:** A GLMM was conducted as in P1-P2.
```{r p5, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
ff <- read.csv ("/Users/corina/GTGR/data/data_flexforaging.csv", header=T, sep=",", stringsAsFactors=F)
#Separate the sexes
fem <- ff[ff$Sex=="f",]
mal <- ff[ff$Sex=="m",]
#Factor the random effect variables
ID <- as.factor(fem$ID)
Population <- as.factor(fem$Population)
ID <- as.factor(mal$ID)
Population <- as.factor(mal$Population)
# DATA CHECKING
library(DHARMa)
library(lme4)
#Data checking for female GLMM
simulationOutput <- simulateResiduals(fittedModel = glmer(ImmigrantProbability ~ TrialsToReverseLast + (1|Population), family=poisson, data=fem), 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(ImmigrantProbability, 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
#Data checking for male GLMM
simulationOutput <- simulateResiduals(fittedModel = glmer(ImmigrantProbability ~ TrialsToReverseLast + (1|Population), family=poisson, data=mal), 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(ImmigrantProbability, 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),R2=list(V=1,nu=0),R3=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0),G2=list(V=1,nu=0)))
#GLMM
#females
f1 <- MCMCglmm(ImmigrantProbability ~ TrialsToReverseLast, random=~Population, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f1)
#autocorr(f1$Sol) #Did fixed effects converge?
#autocorr(f1$VCV) #Did random effects converge?
#males
f2 <- MCMCglmm(ImmigrantProbability ~ TrialsToReverseLast, random=~Population, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f2)
#autocorr(f2$Sol) #Did fixed effects converge?
#autocorr(f2$VCV) #Did random effects converge?
```
####*P6: flexible = wider range of habitats*
**Analysis:** A GLMM was conducted as in P1-P2.
```{r p6, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
ff <- read.csv ("/Users/corina/GTGR/data/data_flexforaging.csv", header=T, sep=",", stringsAsFactors=F)
#Separate the sexes
fem <- ff[ff$Sex=="f",]
mal <- ff[ff$Sex=="m",]
#Factor the random effect variables
ID <- as.factor(fem$ID)
Population <- as.factor(fem$Population)
ID <- as.factor(mal$ID)
Population <- as.factor(mal$Population)
# DATA CHECKING
library(DHARMa)
library(lme4)
#Data checking for female GLMM
simulationOutput <- simulateResiduals(fittedModel = glmer(ShannonDiversityIndexHabitat ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|Population), family=poisson, data=fem), 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(ShannonDiversityIndexHabitat, 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
#Data checking for male GLMM
simulationOutput <- simulateResiduals(fittedModel = glmer(ShannonDiversityIndexHabitat ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition + (1|Population), family=poisson, data=mal), 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(ShannonDiversityIndexHabitat, 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),R2=list(V=1,nu=0),R3=list(V=1,nu=0),R4=list(V=1,nu=0)), G=list(G1=list(V=1,nu=0),G2=list(V=1,nu=0)))
#females
f1 <- MCMCglmm(ShannonDiversityIndexHabitat ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=fem, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f1)
#autocorr(f1$Sol) #Did fixed effects converge?
#autocorr(f1$VCV) #Did random effects converge?
#males
f2 <- MCMCglmm(ShannonDiversityIndexHabitat ~ TrialsToReverseLast + ExperimentalGroup + DominanceRank + Condition, random=~Population+ID, family="poisson", data=mal, verbose=F, prior=prior, nitt=130000, thin=10, burnin=30000)
summary(f2)
#autocorr(f2$Sol) #Did fixed effects converge?
#autocorr(f2$VCV) #Did random effects converge?
```
This species is primarily found within urbanized environments, however there are many different substrates within urban habitats that could provide a variety of food items. Since we are interested in the flexibility of grackle foraging behaviors within the urban habitat, we have focused our habitat diversity measures on the different substrates on which we are mostly likely to see individual variability in foraging behaviors and food types, if present. For example, cement, cafe, and dumpster substrates are all likely to contain human-provided food (either because people leave food out for wild animals or wild animals are able to scrounge human foods), whereas grass, gravel, or other natural substrates such as trees likely contain non-human provided prey items including insects and small vertebrates. Using the Shannon diversity index to understand the evenness of substrate use within urban habitats has been recommended by others in the field of urban ecology (@alberti2001quantifying & @tews2004animal).
###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], SP606267 [2018], SP639866 [2019])
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
**Logan:** Hypothesis development, study design, materials, data collection, data analysis and interpretation, write up, funding.
**Bergeron:** Data collection, data interpretation, revising/editing.
**McCune:** Hypothesis development, study design, data collection, data interpretation, revising/editing.
**Folsom:** Data collection, data interpretation, revising/editing.
**Lukas:** Hypothesis development, study design, data analysis and interpretation, write up, revising/editing.
###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.
###I. 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 Sayres 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; Jay Taylor for grackle scouting at Arizona State University; 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 and helping develop the flexibility 4 estimate; Erin Vogel, our Recommender at PCI Ecology, and Simon Gingins and two anonymous reviewers for their wonderful feedback; two anonymous reviewers and the Recommender, Emanuel Fronhofer, at PCI Ecology for their useful feedback; 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, and Amanda Overholt.
###J. REFERENCES
You can’t perform that action at this time.