Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
2 contributors

Users who have contributed to this file

@corinalogan @kelseybmccune
966 lines (676 sloc) 50.8 KB
---
title: Does great-tailed grackle space use behavior reflect individual differences in exploration?
author: 'McCune KB (University of California Santa Barbara, UCSB, kelseybmccune@gmail.com), Ross C (MPI EVA), Folsom M (Max Planck Institute for Evolutionary Anthropology, MPI EVA), Bergeron L (UCSB), [Logan CJ](http://CorinaLogan.com) (MPI EVA, corina_logan@eva.mpg.de)'
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: /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
Great-tailed grackles (*Quiscalus mexicanus*) are rapidly expanding their geographic range (@wehtje2003range) and it is generally thought that they must rely on behavioral flexibility to achieve this feat [Logan et al. 2019](http://corinalogan.com/Preregistrations/g_flexmanip.html). However, it is alternatively possible that the individuals on the range edge are more exploratory and exhibit distinct movement behaviors (e.g., have larger home ranges and are less predictable about which locations they visit daily). This species is strongly associated with human-modified landscapes and it eats a variety of human foods (e.g., they eat our crops, at our outdoor cafes, and out of our garbage cans) in addition to foraging on insects and in the ground for natural food items (@johnson2001great). The latter hypothesis might be more plausible if the ability to encounter novel foods and food sources is the limiting factor in their expansion rather than their ability to flexibly choose among a variety of options (e.g., to keep track of which restaurants serve lunch at outdoor cafes, when the busiest times are and potentially choose among them according to preferred food type). We aim to understand whether experimental measures of exploration are associated with space use behavior in wild great-tailed grackles in three populations across their range: Central America (their original range), Arizona (middle of the northern expanding edge), and northern California (northern edge of their range). Exploration, measured in the individuals in this study by [McCune et al. 2019](http://corinalogan.com/Preregistrations/g_exploration.html), is interpreted as an individual’s response to novelty, such as novel environments or novel objects (@reale2007integrating), to gather information that does not satisfy immediate needs (@mettke2002significance). We will measure the space use behavior of wild grackles using radio telemetry to find color-banded grackles and record spatial locations across time using GPS. These results will inform whether individual differences in space use behavior are associated with consistent individual differences in exploration, which could be subject to selection, and thus relate to dispersal behaviors (the movement of young and/or adults into new territories) influencing this species’ range expansion within populations. Furthermore, if space use behavior correlates with experimental measures of exploration, then space use could be used to inform conservation management strategies (e.g., which individuals are likely to remain in new or restored habitat after a translocation (@may2016predicting) in species where it is not logistically feasible to experimentally measure exploration. Traditional studies of animal movement behavior require spatial and temporal independence of data points for statistical analysis (@swihart1985testing). However, spatial and temporal autocorrelation (where individuals are found in the same locations across time, such that subsequent relocations are predictable based on previous space use) is an intrinsic component of animal behavior and eliminating it can reduce biological relevance (@dray2010exploratory, e.g., animal movement behavior is influenced by the available habitat and resources which are distributed non-randomly across the landscape). Therefore, in addition to using a typical measure of space use, home range size, we propose two new methods for analyzing wild grackle exploratory behaviors. The first will describe individual differences in movement behavior by analyzing autocorrelation of step length (distance between two sequential observations) and turning angle for each individual over time (@pacheco2019nahua), while the second will describe individual differences in spatial preferences by analyzing the repeatability of each individual's occurrence in particular geographic locations.
###A. STATE OF THE DATA
This preregistration uses secondary data: data that are already being collected for other purposes (GPS points in hypothesis 3 and home range sizes in prediction 3 in the [flexibility and foraging](http://corinalogan.com/Preregistrations/g_flexforaging.html) preregistration). This preregistration was written in June 2019, while at the same time increasing the number of GPS points taken per time per bird to provide enough data for the analyses here, and submitted in September 2019 to PCI Ecology for pre-study peer review.
###B. HYPOTHESES
#### H1: Individual differences in measures of exploration using novel environment and novel object tasks (see separate [preregistration](http://corinalogan.com/Preregistrations/g_exploration.html)) are related to variation in space use (measured via home range size, autocorrelation of step lengths and turning angles, or whether individuals are predictably found in the same locations) across the breeding and non-breeding seasons. During the non-breeding season, this species forages in smaller groups and communally roosts in larger groups [@johnson2001great]. During the breeding season, one or more males defend a territory and females place their nests within territories to raise the young [johnson2000male]. Roaming males are also present and can obtain extra-pair copulations with females on other male's territories [johnson2000male].
**Prediction 1:** The more exploratory grackles that get closer or make more touches to the novel object and novel environment will be found in a larger expanse (home range size), use less predictable movement patterns (low autocorrelation of step lengths), and occupy a greater variety of spatial locations. This suggests that exploratory individuals are more willing to move into novel areas in the wild.
**Prediction 1 alternative 1:** The more exploratory grackles that get closer or make more touches to the novel object and novel environment will use a smaller amount of space (home range size), use predictable movement patterns (high autocorrelation of step lengths), and consistently occupy the same spatial locations. This suggests that more exploratory individuals may be able to more efficiently use habitat within their home range. For example, in great tits, the slow-exploring phenotype relates to more in-depth investigatory behaviors towards changes in the local environment [@verbeek1994consistent], switching to utilization of different resources in the same area [@van2009personality], and better problem-solving abilities [@cole2011personality]. Therefore it may not be necessary for these individuals to move into new areas for resources such as food or mating opportunities.
**Prediction 1 alternative 2:** Only performance on the novel environment task will correlate positively with space use behavior in the wild. This suggests that perception of, and behavioral interactions with, novel environments (spatial information) differs from that used for novel objects [@mettke2009spatial].
**Prediction 1 alternative 3:** Only performance on the novel object task will correlate positively with space use behavior in the wild. This suggests that, in these urban populations, space use may primarily be driven by grackles searching for novel objects that represent human-provided sources of food. Much of the food grackles consume is packaged in human-made packaging (e.g., grackles search in take out bags from restaurants) or enclosed in human-made containers (e.g., garbage cans), therefore they should have a reason to approach and explore new objects to determine whether they could be a new food source.
**Prediction 1 alternative 4:** There will be no correlation between an individual’s proximity or touches to the novel object or novel environment and their space use behavior. This suggests that the experimental measures of exploration either are not relevant enough to how grackles use space in the wild to be able to measure the same trait, or they are independent of space use behavior potentially because the individuals tracked are adults and are already familiar with their home range and surrounding areas and thus do not need to further use the space as if it was novel. Note that only adults are tracked because we only put radio tags on the individuals that are released from the aviaries. We aim to bring only adults in to the aviaries for a cognitive test battery (that is unrelated to this preregistration) so that we are able to understand what this species is capable of, rather than testing juveniles who might still be developing their cognitive skills.
#### H2: Space use behavior will vary among grackles from our three study populations located along different points in the geographic range of this species (core, middle of expansion, and range edge; Table 1).
**Prediction 2:** Home range sizes will increase, autocorrelation of step lengths and turning angles will decrease (grackle movement behavior will be less predictable), and grackles will use a greater variety of spatial locations as the geographic distance from the original center of the range increases. Specifically, grackles on the edge of the range (Northern California), will have larger overall home range sizes, exhibit more variety in step lengths and turning angles, and use a greater variety of spatial locations than grackles in the core of the range (Central America). Grackles in the middle of the expanded range will be intermediate in space use (Arizona). This suggests that range expansion likely occurs as a result of individual differences in space use: the individuals on the leading edge of the range expansion may move longer distances [@duckworth2007coupling].
**Prediction 2 alternative 1:** Grackles on the edge of the range will have smaller overall home range sizes, high autocorrelation in step length and turning angle (movement behavior will be more predictable), and consistently use the same spatial locations compared to grackles in the middle or core of the current range. This suggests that suitable habitat may be distributed in small patches, there may be higher predation on grackles that use more space, and/or individual grackles specialize on certain novel habitat types that are patchily distributed.
**Prediction 2 alternative 2:** There is no difference across the geographic range in the space use behavior of grackles. This suggests that, on average, all grackles use the same amount of space, or that there is a similar distribution of individual differences in space use in each population.
Table 1. Population characteristics for each of the three field sites. Generation length = 5.6 years as estimated by [@birdlife2018Quiscalus].
```{r table1, eval=TRUE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
d <- read.csv ("/Users/corina/GitHub/grackles/Files/Preregistrations/gspaceuse_Table1.csv", header=T, sep=",", stringsAsFactors=F)
library(reactable)
reactable(d, highlight=TRUE, bordered=FALSE, compact=TRUE, wrap=TRUE, resizable=TRUE)
```
###C. METHODS
####**Planned Sample**
Great-tailed grackles are caught in the wild, given colored leg bands in unique combinations for individual identification, and released at their point of capture. The great-tailed grackles used in this study have one of two different backgrounds: those that do not have radio tags and those that do. Those that do have radio tags (estimated 20 individuals per population) are primarily adults who we attempt to balance for sex, and who spent up to 6 months in an aviary while they participated in behavioral choice tests [Logan et al. 2019](http://corinalogan.com/Preregistrations/g_flexmanip.html) and individual differences assays, including captive measures of exploration [McCune et al. 2019](http://corinalogan.com/Preregistrations/g_exploration.html), as part of other research projects by this lab. The radio tags were originally applied to all aviary-tested birds to ensure that we could find their nest sites and track measures of reproductive success for these individuals for which we have extensive amounts of data. This space use preregistration was later developed to address additional questions based on data we were already collecting. For details about the captive environment, please refer to the preregistration associated with this research [McCune et al. 2019](http://corinalogan.com/Preregistrations/g_exploration.html). Before the aviary grackles are released, they are fitted with VHF radio tags (model A2455 by Advanced Telemetry Systems or model BD-2 by Holohil Systems Ltd.) so we can track space use behavior using radio telemetry. Radio tags were initially attached to the grackles by gluing them to their backs (@johnson2001great, @mong2007optimizing), however these did not stay on for very long. Therefore, now we use a leg loop harness (methods as in @rappole1991new) made from sutures (Vicryl undyed 36in sutures, item number D9389 at eSutures.com; 0.5mm diameter, absorbable so they fall off after a few months).
After release, an experimenter tracks each tagged grackle for approximately 1.5 hours on a given day, taking a GPS point every 1 minute or every 10 minutes (we will determine the sampling interval that gives the best resolution), regardless of whether the bird moved [@cushman2005elephants]. Researchers are careful not to get too close so as not to influence the grackle’s behavior.
####**Sample size rationale**
We test as many birds as we can in the approximately five years of this study given that the birds are only brought into the aviaries during the non-breeding season (approximately September through March). It is time intensive to conduct the aviary test battery (2-6 months per bird at the Arizona field site), therefore we approximate that the minimum sample size for captive subjects will be 57 across the three sites (approximately 20 birds per site with half of the grackles at each site being female). 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). Once released, we primarily track the space use behavior of the grackles that were in the aviary, but we also collect GPS point locations on all occasions that we see any color-marked bird so we can determine whether grackles that were previously in the aviary have different space use behavior from non-aviary-held grackles after their release.
####**Data collection stopping rule**
We will stop collecting GPS location data on tagged birds when home ranges are fully revealed. To determine at what point home ranges have been fully revealed, we will calculate the asymptotic convergence of home range area as in @leo2016home. We will test home range asymptotic convergence for breeding season and non-breeding season movements separately (breeding season: Apr - Aug, non-breeding season: Sep - Mar).
#### **Open materials**
Testing protocols:
- [Exploration protocol](https://docs.google.com/document/d/1sEMc5z2fw6S9C-wVfc2zV331CRPpu3NuA7IhSFUZJpE/edit?usp=sharing) for exploration of new environments and objects, boldness, persistence, and motor diversity.
- [Radio tracking protocol](https://docs.google.com/document/d/1jtjgeWJoZ0Q1CfUpV6zdkyQL3p3WfW9KgyLrMNmNMJc/edit?usp=sharing) for attaching radio tags and collecting GPS points using radio telemetry.
####**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 randomization in this investigation. The order of the exploration tasks is counterbalanced across birds as in this [separate preregistration](http://corinalogan.com/Preregistrations/g_exploration.html). The time of day that we collect GPS point locations is counterbalanced within and across birds to account for variation in movement behavior arising from daily circadian rhythms.
####**Blinding of conditions during analysis**
No blinding is involved in this investigation.
####**Summary of methods for measuring exploration [McCune et al. 2019](http://corinalogan.com/Preregistrations/g_exploration.html)**
Exploration assays for the grackles that are temporarily held in aviaries occur twice: once near the beginning of their aviary time and once again approximately 6 weeks later. We will analyze whether behavioral responses during assays are repeatable within individuals. The order of assays (exploration novel environment, exploration novel object, and boldness) is counterbalanced across individuals.
####**Dependent variables**
**P1-P2**
1) Home range size (square meters): an estimate calculated using the autocorrelated-Gaussian reference function kernel density estimate (AKDE), which is the only estimate of home range that accounts for the autocorrelation due to the small time period between each or our GPS locations [@noonan2019comprehensive]. This estimate consists of the area enclosing 95% of GPS location points for an individual grackle during its normal activities.
2) Autocorrelation of step length (meters): measured as the standard deviation of step lengths (the distance between two sequential GPS points)
3) Autocorrelation of turning angle (degrees): measured as the standard deviation of turning angles
4) Spatial location preference: measured as the repeatability of grackle occurrence in a given cell of a 5 x 5m grid array across the landscape
One model will be run for each dependent variable
####**Independent variables**
**P1 and P1 alternatives 2-4**
1) Exploration of novel environment: Latency to approach to 20cm of a novel environment (that does not contain food) set inside a familiar environment (that contains maintenance diet away from the object) - OR - closest approach distance to the novel environment (choose the variable with the most data)
2) Exploration of novel object: Latency to approach to 20cm of an object (novel or familiar, that does not contain food) in a familiar environment (that contains maintenance diet away from the object) - OR - closest approach distance to the object (choose the variable with the most data)
3) Sex: Male or female
4) History: the number of days the individual was temporarily held in the aviaries before data collection on space use began (0 indicates the grackle was only ever in the wild)
**P1 alternative 1**
5) Number of different foods eaten and foraging techniques used: Measured during focal follow observations on wild grackles, described in this [separate preregistration](http://corinalogan.com/Preregistrations/g_flexforaging.html)
**P2**
1) Site: Whether the grackle was from our study population located on the edge of the range (Northern California), the center of the original range (Central America), or the center of the current expanded range (Arizona).
2) Sex: Male or female
3) History: the number of days the individual was temporarily held in the aviaries before data collection on space use began (0 indicates the grackle was only ever in the wild)
###D. ANALYSIS PLAN
We do not plan to **exclude** any data and if there are **missing** data (e.g., if a bird had to be released before collecting their data at time 2) these birds will be excluded from analyes requiring data from time 1 and 2. Analyses will be conducted in R (current version `r getRversion()`; @rcoreteam) and Stan (version 2.18, @carpenter2017stan).
We will first verify that the GPS point locations on each bird result in asymptotic convergence as in @leo2016home. To calculate our dependent variables we will use the autocorrelated kernel density estimate method for quantifying home range size (in square meters) using the akde function in the R packages ctmm [@calabrese2016ctmm] and sf [@pebesma2018simple]. Autocorrelated kernel density estimates (AKDE) of home range size are the most accurate when data are collected close together in time and space [@noonan2019comprehensive]. We are interested in all movements by grackles, so we will not exclude any outlier relocations collected during “normal activities” [@calenge2011home]. ‘Normal activities’ indicate that grackles are not engaging in behaviors that would artificially skew their space use, for example mobbing a predator or the experimenter, or behavior before sunrise or after sunset when they are at the roost. Outside of these circumstances, we will include all data to detect any space use movements.
From the GPS point locations collected on each individual, we will use a Bayesian model (detailed below) to estimate the following parameters: mean and dispersion (variance) of step lengths and turning angles for each bird on each daily track (@pacheco2019nahua). We will determine whether these parameters governing movement are stable (or variable) within individuals across days. A small variance would indicate there is low variability (high repeatability) in the daily movement behaviors of the individual.
Moreover, we will determine whether grackles show individual differences in consistent use of habitat by overlaying a grid array across the landscape. We will then create matrices describing the number of times a grackle was observed in each cell on each day. High autocorrelation among daily matrices indicates a lesser moving individual that frequents the same spatial locations across days.
We will then model the relationship between bird-specific data on performance in our exploration tasks (and other covariates), and bird-specific movement parameters (e.g. step-size, turning angle, autocorrelation in space use).
####*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). We realize that these power analyses are not fully aligned with our study design, 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.
#### *Calculating home range size*
```{r hr, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#Load packages
library(ctmm)
library(sf)
#Point to the correct data file and load it
setwd("~/Documents/Grackle project/Space use")
pts<-read.csv("gtgr_points.csv", header = T)
#Set variables to ensure the models read the data properly
pts$Latitude = as.numeric(as.character(pts$Latitude))
pts$Longitude = as.numeric(as.character(pts$Longitude))
pts$Date = as.character(pts$Date)
pts$Time = as.character(pts$Time)
pts$da.ti = as.POSIXct(paste(pts$Date, pts$Time), format="%Y-%m-%d %H:%M:%S")
### Home range size can only be calculated when a bird has more than 5 relocations
### Count the number of relocations for each bird, then exclude individuals with less than 6
tmp = pts
tmp$count = 1
tmp = aggregate(count ~ Bird.Name, FUN = "sum", data = tmp)
pts.min = merge(pts, tmp, by = "Bird.Name", all = T)
pts.min = pts.min[-which(pts.min$count<6),]
pts.min$Bird.Name = as.factor(as.character(pts.min$Bird.Name))
### Convert dataframe to Spatial file type
p.sf <- st_as_sf(pts.min, coords = c("Longitude", "Latitude"), crs = 4326)
class(p.sf)
p.spatial <- as(p.sf, "Spatial")
class(p.spatial)
### Calculate home range in square meters for each individual
hr = mcp(p.spatial[1], percent = 100, unout = "m2")
plot(hr)
hr = as.data.frame(hr)
### Normalize data for regression model
hist(log(hr$area))
### Make data sheet with the experimental exploration measures combined with the home range measures
exp <- read.csv("Explore_combined.csv", header = T)
hr_exp = merge(exp, hr, by = "id", all = T)
write.csv(hr_exp, "space_use.csv")
```
#### *Code to create functions for analyzing movement behaviors*
All scripts and code are available at https://github.com/ctross/grackleator.
```{r functions, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
### gracklenomial.R
gracklenomial <- function(GrackBins){
model_dat_2=list(
Trips=dim(GrackBins)[1],
Bins=dim(GrackBins)[2],
GrackBins=GrackBins)
model_code_2 <- '
data{
int Trips;
int Bins;
int GrackBins [Trips,Bins];
}
parameters{
simplex[Bins] A;
real<lower=0,upper=1> B;
}
model{
vector[Bins] C;
A ~ dirichlet(rep_vector(1,Bins));
B ~ beta(1,1);
for( i in 2:Trips){
C = to_vector(GrackBins[i-1])/sum(GrackBins[i-1]);
GrackBins[i] ~ multinomial( A*(1-B) + B*C );
}
}
'
m2 <- stan( model_code=model_code_2, data=model_dat_2,refresh=1,chains=1, control = list(adapt_delta = 0.9, max_treedepth = 13))
return(m2)
}
### gracklepars.R
gracklepars <- function(m1){
MuD_M<-data.frame(Group="Step-Size",Group2="Mean",Group3="Mean",Value=rstan::extract(m1,pars="MuD_M")$MuD_M) # Mean Step-Size over Days
DispD_M<-data.frame(Group="Step-Size",Group2="Mean",Group3="Dispersion",Value=rstan::extract(m1,pars="DispD_M")$DispD_M) # Dispersion in Mean Step-Size over days
MuD_D <-data.frame(Group="Step-Size",Group2="Dispersion",Group3="Mean",Value=rstan::extract(m1,pars="MuD_D")$MuD_D)# Mean Dispersion in Step-Size
DispD_D <-data.frame(Group="Step-Size",Group2="Dispersion",Group3="Dispersion",Value=rstan::extract(m1,pars="DispD_D")$DispD_D)# Dispersion in Dispersion in Step-Size
MuA_M <-data.frame(Group="Heading Change",Group2="Mean",Group3="Mean",Value=rstan::extract(m1,pars="MuA_M")$MuA_M)# Mean Angle Change over Days
DispA_M <-data.frame(Group="Heading Change",Group2="Mean",Group3="Dispersion",Value=rstan::extract(m1,pars="DispA_M")$DispA_M)# Dispersion in Mean Angle Change over days
MuA_D <-data.frame(Group="Heading Change",Group2="Dispersion",Group3="Mean",Value=rstan::extract(m1,pars="MuA_D")$MuA_D)# Mean Dispersion in Angle Change
DispA_D <-data.frame(Group="Heading Change",Group2="Dispersion",Group3="Dispersion",Value=rstan::extract(m1,pars="DispA_D")$DispA_D)# Dispersion in Dispersion in Angle Change
df_allx<-rbind(MuD_M,DispD_M,MuD_D,DispD_D,MuA_M,DispA_M,MuA_D,DispA_D)
ggplot(df_allx, aes(x=Value)) + geom_density(aes(group=Group3, colour=Group3, fill=Group3), alpha=0.3) + facet_wrap(Group2~Group,scales="free")+
labs(y="Regression parameters") + theme(strip.text.x = element_text(size=14,face="bold"),axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))
}
### grackletrip.R
grackletrip <- function(m1){
m1a<-rstan::extract(m1,pars="AlphaAngle")
sample_eff<-apply(m1a$AlphaAngle,2,quantile,probs=c(0.05,0.5,0.95))
df_angle<-data.frame(Trip=c(1:dim(m1a$AlphaAngle)[2]),Group="Heading Change",Group2="Mean",
LI=sample_eff[1,],
Median=sample_eff[2,],
HI=sample_eff[3,])
m1d<-rstan::extract(m1,pars="AlphaDist")
sample_eff<-apply(m1d$AlphaDist,2,quantile,probs=c(0.05,0.5,0.95))
df_dist<-data.frame(Trip=c(1:dim(m1d$AlphaDist)[2]),Group="Step-Size",Group2="Mean",
LI=sample_eff[1,],
Median=sample_eff[2,],
HI=sample_eff[3,])
df_all1<-rbind(df_angle,df_dist)
m1a<-rstan::extract(m1,pars="DAngle")
sample_eff<-apply(exp(m1a$DAngle),2,quantile,probs=c(0.05,0.5,0.95))
df_angle<-data.frame(Trip=c(1:dim(m1a$DAngle)[2]),Group="Heading Change",Group2="Dispersion",
LI=sample_eff[1,],
Median=sample_eff[2,],
HI=sample_eff[3,])
m1d<-rstan::extract(m1,pars="SDDist")
sample_eff<-apply(exp(m1d$SDDist),2,quantile,probs=c(0.05,0.5,0.95))
df_dist<-data.frame(Trip=c(1:dim(m1d$SDDist)[2]),Group="Step-Size",Group2="Dispersion",
LI=sample_eff[1,],
Median=sample_eff[2,],
HI=sample_eff[3,])
df_all2<-rbind(df_angle,df_dist)
df_all<-rbind(df_all1,df_all2)
ggplot(df_all,aes(x=Trip,y=Median))+geom_point()+
geom_linerange(aes(ymin=LI,ymax=HI))+facet_wrap(Group2~Group,scales="free")+
labs(y="Regression parameters") + theme(strip.text.x = element_text(size=14,face="bold"),axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold")) +
scale_x_continuous(breaks= pretty_breaks())
}
### grackleate.R
grackleate <- function(AlphaDist=2.8,AlphaAngle=0,BetaDist=rep(0,10),BetaAngle=rep(0,10),SDDist=1,DAngle=2,Thresh=50, Lags=10, steps=150, seed=123456, Reps=30,N_Patches=13){
################################################################################# Set up Patches of Food
X_Mushrooms<-vector("list",N_Patches) # Location of prey
Y_Mushrooms<-vector("list",N_Patches) #
set.seed(seed) # Reset Seed
N_PerPatch<-rpois(N_Patches, 200) # Create some prey
for(i in 1:N_Patches){
X_Mushrooms[[i]]<-rnorm(N_PerPatch[i], runif(1,-2500,2500), rpois(1, 350))+100
Y_Mushrooms[[i]]<-rnorm(N_PerPatch[i], runif(1,-2500,2500), rpois(1, 350))+100
}
X_Mushrooms_per_step<-vector("list",steps) # Location of prey
Y_Mushrooms_per_step<-vector("list",steps) #
for(i in 1:steps){
X_Mushrooms_per_step[[i]]<-X_Mushrooms
Y_Mushrooms_per_step[[i]]<-Y_Mushrooms
}
################################################################################ Start Model
Store.X <- matrix(NA,ncol=Reps,nrow=steps-1)
Store.Y <- matrix(NA,ncol=Reps,nrow=steps-1)
for(q in 1:Reps){ # Run simulation several times
set.seed(seed) # Reset seed
N_PerPatch<-rpois(N_Patches, 200) # Choose number of items per patch
for(i in 1:N_Patches){
X_Mushrooms[[i]]<-rnorm(N_PerPatch[i], runif(1,-2500,2500), rpois(1, 350))+100 # Make some patches of prey
Y_Mushrooms[[i]]<-rnorm(N_PerPatch[i], runif(1,-2500,2500), rpois(1, 350))+100 #
}
X_Mushrooms_per_step<-vector("list",steps) # Location of prey
Y_Mushrooms_per_step<-vector("list",steps) #
for(i in 1:steps){
X_Mushrooms_per_step[[i]]<-X_Mushrooms
Y_Mushrooms_per_step[[i]]<-Y_Mushrooms
}
loc.x<-c() # Location of forager
loc.y<-c() #
speed<-c() # Speed or step size
heading<-c() # Absolute heading
d.heading<-c() # Heading change
Hits<-c() # Binary vector of encounters
loc.x[1:Lags]<-rep(0,Lags) # Intialize vectors
loc.y[1:Lags]<-rep(0,Lags) #
heading[1:Lags]<-rep(0,Lags) #
speed[1:Lags]<-rep(0,Lags) #
Hits[1:Lags]<-rep(0,Lags) #
plot(loc.x,loc.y,typ="l",ylim=c(-3500,3500),xlim=c(-3500,3500)) # Plot prey items
for(i in 1:N_Patches){ #
points(X_Mushrooms[[i]],Y_Mushrooms[[i]], col="red",pch=".") #
} #
set.seed(q*103)
################################################################################ Now model forager movement
for(s in (Lags+1): (steps-1)){
X_Mushrooms <- X_Mushrooms_per_step[[s]]
Y_Mushrooms <- Y_Mushrooms_per_step[[s]]
PredDist <- AlphaDist; # First calculate mean prediction conditional on encounters
for(k in 1:Lags){ #
PredDist <- PredDist + BetaDist[k]*ifelse(Hits[s-k]>0,1,0); #
} #
R<- exp(rnorm(1,PredDist,SDDist)) # Then simulate a step distance.
PredAngle <- AlphaAngle; # Again calculate mean prediction conditional on encounters
for(k in 1:Lags){ #
PredAngle <- PredAngle + BetaAngle[k]*ifelse(Hits[s-k]>0,1,0); #
} #
Theta<- rbeta(1,inv_logit(PredAngle)*DAngle, # And then simulate a directional change
(1-inv_logit(PredAngle))*DAngle)*180*ifelse(runif(1,0,1)>0.5,1,-1) # The 180 shifts from the unit to the maximum absolute distance
# The ifelse just chooses a random direction
heading[s]<-(heading[s-1]+Theta)%%360 # Store new heading. Note that the %% is the mod operation to wrap around if needed
d.heading[s] <- abs(Theta)/180 # Also store just the delta heading
speed[s] <- R # And the speed slash step size
ynew <- R * sin(deg2rad(heading[s])) # Now convert polar to Cartesian, to get the offset for a new x and y pair
xnew <- R * cos(deg2rad(heading[s])) #
loc.x[s]<-loc.x[s-1]+xnew # Make the new x and y pair
loc.y[s]<-loc.y[s-1]+ynew #
############################################################################### Now check for an encounter
Scrap2<-c()
for(i in 1:N_Patches){ # For each patch
Scrap<-rep(NA,length(X_Mushrooms[[i]]))
for(j in 1:length(X_Mushrooms[[i]])){ # For each mushroom in patch
Scrap[j]<- dist(loc.x[s],X_Mushrooms[[i]][j],loc.y[s],Y_Mushrooms[[i]][j]); # Calculate the distance from the forager to the mushroom
if(Scrap[j]<Thresh){ # If the forager is closer than the visual radius slash encounter threshold
points(X_Mushrooms[[i]][j],Y_Mushrooms[[i]][j], col="blue",pch=20) # Plot a hit
for(allgone in s: (steps-1)){
X_Mushrooms_per_step[[allgone]][[i]][j]<-99999 # If this run is for destructive foraging
X_Mushrooms_per_step[[allgone]][[i]][j]<-99999 # disappear food for all future timesteps
}
}}
Scrap2[i]<- ifelse(sum(ifelse(Scrap<Thresh,1,0))>0,sum(ifelse(Scrap<Thresh,1,0)),0) # Check for hits, also can replace sum(ifelse(Scrap<Thresh,1,0)) with 1
}
Hits[s]<-ifelse(sum(Scrap2)==0,0,sum(Scrap2)) # If there is a hit, then set encounters to 1
lines(loc.x,loc.y,ylim=c(-2500,2500),xlim=c(-2500,2500)) # Plot updates to the foragers path
}
Store.X[,q] <- loc.x
Store.Y[,q] <- loc.y
}
return(list(X=Store.X[(Lags+1):length(loc.x),],Y=Store.Y[(Lags+1):length(loc.y),]))
}
### grackleations.R
grackleations <- function(m2){
df2<-data.frame(Parameter="B",Value=rstan::extract(m2,pars="B")$B) # Mean Step-Size over Days
ggplot(df2, aes(x=Value)) + geom_density(aes(fill=Parameter), alpha=0.3) +
theme(strip.text.x = element_text(size=14,face="bold"),axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))
}
### gracklebinner.R
gracklebinner <- function(tracks,nbin=c(15,15),ab_override=NULL){
Trips <- dim(tracks$X)[2]
GrackBins <- matrix(NA,nrow=Trips,ncol=nbin[1]*nbin[2])
bins<-bin2(cbind(c(z$X),c(z$Y)),nbin=nbin)
if(length(dim(ab_override))==0){
for( i in 1:Trips){
GrackBins[i,] <- c(bin2(cbind(z$X[,i],z$Y[,i]),nbin=nbin,ab=bins$ab)$nc)
}
} else{
for( i in 1:Trips){
GrackBins[i,] <- c(bin2(cbind(z$X[,i],z$Y[,i]),nbin=nbin,ab=ab_override)$nc)
}
}
return(GrackBins)
}
### grackleize.R
grackleize <- function(StoreX,StoreY) {
############################################################ Prepare Data
Trip <- dim(StoreX)[2]
MaxTicks<-dim(StoreX)[1]
Distance <-matrix(NA,ncol=max(Trip),nrow=MaxTicks)
Ang <-matrix(NA,ncol=max(Trip),nrow=MaxTicks)
AngDiff <-matrix(NA,ncol=max(Trip),nrow=MaxTicks)
N<-c()
for(j in 1:max(Trip)){
X <- StoreX[,j]
Y <- StoreY[,j]
N[j] <- length(X)
Dist<-c()
ang1<-c()
angdif<-c()
for(i in 2:N[j]){
Dist[i]<- dist(X[i],X[i-1],Y[i],Y[i-1]);
ang1[i]<- ang(X[i],X[i-1],Y[i],Y[i-1])
if(i>3){
angdif[i]<- ang.dif(ang1[i],ang1[i-1])
}}
Distance[1:N[j],j]<-Dist
Ang[1:N[j],j]<-ang1
AngDiff[1:N[j],j]<-angdif
}
for(i in 1:MaxTicks){
for(j in 1:max(Trip)){
Distance[i,j] <- ifelse(Distance[i,j]==0,
runif(1,0,min(Distance[which(Distance != 0)],na.rm=TRUE)),
Distance[i,j])
}}
Distance[is.na(Distance)] <- 999999
EmptyAngle <- fitdistr( ScaleBeta(c(na.omit(c(AngDiff)))),start=list(shape1=1,shape2=1),"beta")$estimate
for(i in 1:MaxTicks){
for(j in 1:max(Trip)){
AngDiff[i,j] <- ifelse(is.nan(AngDiff[i,j]), rbeta(1,EmptyAngle[1], EmptyAngle[2])*pi,AngDiff[i,j])
}}
AngDiff <- ScaleBeta(AngDiff)
AngDiff[is.na(AngDiff)] <- 999999
AngDiff <- AngDiff[c(-1,-2,-3),]
Distance <- Distance[c(-1,-2,-3),]
N <- N-3
MaxTicks<-MaxTicks-3
########################################################## Extract data for Stan
model_dat=list(
Distance=Distance,
AngDiff=AngDiff,
N=N,
MaxTrip=max(Trip),
MaxTicks=MaxTicks)
################################################################# Fit Stan Model
model_code_1 <- '
data{
int MaxTrip;
int MaxTicks;
int N[MaxTrip];
real Distance[MaxTicks,MaxTrip];
real AngDiff[MaxTicks,MaxTrip];
}
parameters {
real MuD_M;
real MuA_M;
real<lower=0> DispD_M;
real<lower=0> DispA_M;
real MuD_D;
real MuA_D;
real<lower=0> DispD_D;
real<lower=0> DispA_D;
vector[MaxTrip] AlphaDist_Raw;
vector[MaxTrip] SDDist_Raw;
vector[MaxTrip] AlphaAngle_Raw;
vector[MaxTrip] DAngle_Raw;
}
transformed parameters{
vector[MaxTrip] AlphaDist;
vector[MaxTrip] SDDist;
vector[MaxTrip] AlphaAngle;
vector[MaxTrip] DAngle;
AlphaDist = MuD_M + DispD_M*AlphaDist_Raw;
SDDist = MuD_D + DispD_D*SDDist_Raw;
AlphaAngle = MuA_M + DispA_M*AlphaAngle_Raw;
DAngle = MuA_D + DispA_D*DAngle_Raw;
}
model{
MuD_M~normal(0,5);
MuA_M~normal(0,5);
DispD_M~cauchy(0,1);
DispA_M~cauchy(0,1);
MuD_D~normal(0,5);
MuA_D~normal(0,5);
DispD_D~cauchy(0,1);
DispA_D~cauchy(0,1);
AlphaDist_Raw~normal(0,1);
SDDist_Raw~normal(0,1);
AlphaAngle_Raw~normal(0,1);
DAngle_Raw~normal(0,1);
for(j in 1:MaxTrip){
{
vector[N[j]] PredAngle;
vector[N[j]] PredDist;
vector[N[j]] Dist;
vector[N[j]] AngleDifference;
for(i in 1:N[j]){
PredDist[i] = AlphaDist[j];
}
for(i in 1:N[j]){
PredAngle[i] = AlphaAngle[j];
}
for(i in 1:N[j]){
Dist[i] = Distance[i,j];
AngleDifference[i] = AngDiff[i,j];
}
Dist ~ lognormal(PredDist,exp(SDDist[j]));
AngleDifference ~ beta(inv_logit(PredAngle)*exp(DAngle[j]), (1-inv_logit(PredAngle))*exp(DAngle[j]));
}}
}
'
m1 <- stan( model_code=model_code_1, data=model_dat,refresh=1,chains=1, control = list(adapt_delta = 0.9, max_treedepth = 15))
return(m1)
}
### setup.R
# Load libraries
library(MASS)
library(mvtnorm)
library(fields)
library(ggplot2)
library(rethinking)
library(sfsmisc)
library(ash)
library(reshape)
library(maptools)
library(gridExtra)
library(scales)
library(msm)
library(maps)
library(grid)
library(xtable)
# Define functions
dist2<-function(a,b){return(sqrt( (b[2]-a[2])^2 + (b[1]-a[1])^2 ))}
dist<-function(x2a,x1a,y2a,y1a){return(sqrt((x2a-x1a)^2 + (y2a-y1a)^2))}
rad2deg <- function(rad) {(rad * 180) / (pi)}
deg2rad <- function(deg) {(deg * pi) / (180)}
ang.dif <- function(x,y) {min((2 * pi) - abs(x - y), abs(x - y))}
lp_dist <-function(a,b,c){
# Return distance between a point and line segment
t <- -(((a[1]-c[1])*(b[1]-a[1]) + (a[2]-c[2])*(b[2]-a[2])) / ((b[1]-a[1])^2 + (b[2] - a[2])^2 ))
if(t >0 & t <1){
numer <- abs( (b[2]-a[2])*c[1] -(b[1]-a[1])*c[2] + b[1]*a[2] - b[2]*a[1])
denom <- sqrt( (b[2]-a[2])^2 + (b[1]-a[1])^2 )
return(numer/denom)
}
else{
d1 <- dist2(a,c)
d2 <- dist2(b,c)
return( min(d1,d2))
}
}
ScaleBeta <- function(X){
a<-0
b<-pi
y<-X
Samp<-50
y2 <- (y-a)/(b-a)
y3<-(y2*(Samp - 1) + 0.5)/Samp
y3}
ang<-function(x2,x1,y2,y1){
Theta<-seq(0,2*pi,length.out=100)
DB<-cbind(cos(Theta),sin(Theta))
r <- dist(x2,x1,y2,y1)
x0<-(x2-x1)/r
y0<-(y2-y1)/r
(atan2(y0,x0) + pi )
}
################################################################# Set parameters
Thresh <- 50 # Threshold of visual range
Lags<-10 # Lags at which there are effects of encounters on movement, kinda hardcoded---careful with changes
Steps<-150 # Length of simulation
Reps <- 30 # Number of replications
# Adaptive model parameters
Phi0A <- 3.2
Psi0A <- 0
PhiA <- c(-0.76, -0.65, -0.58, -0.43, -0.29, -0.15, -0.03,0,0,0)
PsiA <- c(0.5, 0.3, 0.25, 0.2, 0.17, 0.15, 0.13, 0.11, 0.09, 0.08)
OmegaA <- 0.4
EtaA <- 2
# Levy parameters
Phi0L <- 2.8
Psi0L <- 0
PhiL <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
PsiL <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
OmegaL <- 1
EtaL <- 2
# Brownian parameters
Phi0B <- 25
Psi0B <- 0
PhiB <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
PsiB <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
OmegaB <- 15
EtaB <- 2
```
#### *Modeling bird movement behaviors: step length, turning angle, spatial location preference*
```{r grackleator, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
#This is an R package for modeling bird movement (available at https://github.com/ctross/grackleator)
#Install by running on R:
library(devtools)
install_github("Ctross/grackleator")
#Some quick examples:
#First, lets look at movement dynamics:
# Load library and attach data
library(grackleator)
# Simulate tracks from 1 grackle over 30 trips
z <- grackleate(AlphaDist=3.8,AlphaAngle=0,SDDist=1.5,DAngle=2)
# Analyze results of trips for step-size and angle change
m1 <- grackleize(z$X,z$Y)
# Plot results over trips
grackletrip(m1)
# Plot bird-specific parameters
gracklepars(m1)
#And now habitat selection:
# Load library and attach data
library(grackleator)
# Simulate tracks from 1 grackle over 30 trips
z <- grackleate(AlphaDist=3.8,AlphaAngle=0,SDDist=1.5,DAngle=2)
# Bin data on 2D grid
GrackBins<- gracklebinner(z,ab_override=matrix(c(-3000,-3000,3000,3000),nrow=2,ncol=2))
# Model location data
m2 <- gracklenomial(GrackBins)
# Plot environmental hotspots
image(matrix(colSums(GrackBins), nrow=15,ncol=15)) # Overall
image(matrix(GrackBins[1,], nrow=15,ncol=15)) # Day 1
image(matrix(GrackBins[2,], nrow=15,ncol=15)) # Day 2
image(matrix(GrackBins[3,], nrow=15,ncol=15)) # Day 3
image(matrix(get_posterior_mean(m2,pars="A"), nrow=15,ncol=15)) # Overall
# Plot bird-specific parameters
grackleations(m2)
# Now create data with temporal correlations
GrackBins2 <- GrackBins
A <- get_posterior_mean(m2,pars="A")
A <- A/sum(A)
B <- 0.8
for(i in 2:30)
GrackBins2[i,] <- rmultinom(1,sum(GrackBins2[i-1,]),A*(1-B) + B*(GrackBins2[i-1,]/sum(GrackBins2[i-1,])))
# Analyze the new data
m3 <- gracklenomial(GrackBins2)
# Plot environmental hotspots
image(matrix(colSums(GrackBins2), nrow=15,ncol=15)) # Overall
image(matrix(GrackBins2[1,], nrow=15,ncol=15)) # Day 1
image(matrix(GrackBins2[2,], nrow=15,ncol=15)) # Day 2
image(matrix(GrackBins2[3,], nrow=15,ncol=15)) # Day 3
image(matrix(get_posterior_mean(m3,pars="A"), nrow=15,ncol=15)) # Post
# Plot bird-specific parameters
grackleations(m3)
```
#### *H1: P1 - Exploration measured in captivity relates to space use behavior*
To roughly estimate our ability to detect actual effects, we ran a power analysis in G\*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,19
α err prob = 0,05
Power (1-β err prob) = 0,7
Number of predictors = 4
*Output:*
Noncentrality parameter λ = 10,6400000
Critical F = 2,5533954
Numerator df = 4
Denominator df = 51
Total sample size = 56
Actual power = 0,7009879
This means that, with our minimum sample size of 57, we have a 70% chance of detecting a small effect (approximated at f^2^=0.20 by @cohen1988statistical).
```{r h1, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
data <- read.csv("Space_use.csv", header = T)
# Home range
m1 = lm(log(area) ~ ExpObj + ExpEnv + Sex + Condition, data = data)
hist(m1$resid)
summary(m1)
# Step length
m2 = lm(log(std_step)) ~ ExpObj + ExpEnv + Sex + Condition, data = data)
hist(m2$resid)
summary(m2)
# Turning angle
m3 = lm(log(std_angle)) ~ ExpObj + ExpEnv + Sex + Condition, data = data)
hist(H3$resid)
summary(m3)
# Spatial preferences
m4 = lm(log(loc_pref)) ~ ExpObj + ExpEnv + Sex + Condition, data = data)
hist(H4$resid)
summary(m4)
```
#### *H2: P2 - Space use behaviors vary among populations across the range*
To roughly estimate our ability to detect actual effects, we ran a power analysis in G\*Power in the same way as for Hypothesis 1. The protocol of the power analysis is here:
*Input:*
Effect size f² = 0,17
α err prob = 0,05
Power (1-β err prob) = 0,7
Number of predictors = 3
*Output:*
Noncentrality parameter λ = 9,5200000
Critical F = 2,7826004
Numerator df = 3
Denominator df = 52
Total sample size = 56
Actual power = 0,7021017
This means that, with our minimum sample size of 57, we have a 70% chance of detecting a small effect (approximated at f^2^=0.20 by @cohen1988statistical).
```{r h2, eval=FALSE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
data <- read.csv("Space_use.csv", header = T)
# Home range
m1 = lm(log(area) ~ Site + Sex + Condition, data = data)
hist(m1$resid)
summary(m1)
# Step length
m2 = lm(log(std_step)) ~ Site + Sex + Condition, data = data)
hist(m2$resid)
summary(m2)
# Turning angle
m3 = lm(log(std_angle)) ~ Site + Sex + Condition, data = data)
hist(m3$resid)
summary(m3)
# Spatial preference
m4 = lm(log(loc_pref)) ~ Site + Sex + Condition, data = data)
hist(m4$resid)
summary(m4)
```
###E. 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], and SP639866 [2019])
4) Institutional Animal Care and Use Committee at Arizona State University (protocol number 17-1594R)
###F. AUTHOR CONTRIBUTIONS
**McCune:** Hypothesis development, data collection, data analysis and interpretation, write up, revising/editing.
**Folsom:** Data collection, revising/editing.
**Ross:** Model development, data analysis and interpretation, revising/editing.
**Bergeron:** Data collection, revising/editing.
**Logan:** Hypothesis development, data interpretation, write up, revising/editing, materials/funding.
###G. FUNDING
This research is funded by the Department of Human Behavior, Ecology and Culture at the Max Planck Institute for Evolutionary Anthropology.
###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 Melissa Wilson Sayres for sponsoring our affiliations at Arizona State University (ASU); Kevin Langergraber for serving as local PI on the ASU IACUC; Kristine Johnson for technical advice on great-tailed grackles; Julia Cissewski for tirelessly solving problems involving financial transactions and contracts; Richard McElreath for project support; and our research assistants: Aldora Messinger, Elysia Mamola, Michael Guillen, Rita Barakat, Adriana Boderash, Olateju Ojekunle, August Sevchik, Justin Huynh, Jennifer Berens, Michael Pickett and Amanda Overholt.
###K [REFERENCES](MyLibrary.bib)
You can’t perform that action at this time.