gspaceuse.Rmd

---
title: Does space use behavior relate to exploration in a species that is rapidly expanding its geographic range?
author: 
- '[McCune KB](https://www.kelseymccune.com/)^1^*'
- 'Ross C^2^'
- 'Folsom M^2^'
- 'Bergeron L^1^'
- '[Logan CJ](http://CorinaLogan.com)^2^'
date: '`r Sys.Date()`'
output:
  html_document: 
    toc: true
    toc_depth: 4
    toc_float: 
      collapsed: false
    code_folding: hide 
  md_document: 
    toc: true
  pdf_document:
    keep_tex: yes
    latex_engine: xelatex
  github_document: 
    toc: true
bibliography: /Users/corina/ownCloud/Documents/GitHub/grackles/Files/MyLibrary.bib
---

##### Affiliations: 
1) University of California Santa Barbara
2) Max Planck Institute for Evolutionary Anthropology

*Corresponding author: KB McCune (kelseybmccune@gmail.com)

```{r setup, include=FALSE}
library(knitr)
knitr::opts_chunk$set(echo = TRUE)

#Make code wrap text so it doesn't go off the page when Knitting to PDF
opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=TRUE)

#Trying to get Cody's code to knit but Tidy is messing it up https://github.com/ramnathv/poirot/issues/25
opts_chunk$set(tidy = F)
```

<img width="50%" src="logoPCIecology.png">

**Cite as:** McCune K, Ross C, Folsom M, Bergeron L, Logan CJ. 2020. [Does space use behavior relate to exploration in a species that is rapidly expanding its geographic range?](http://corinalogan.com/Preregistrations/gspaceuse.html) (http://corinalogan.com/Preregistrations/gspaceuse.html) In principle acceptance by *PCI Ecology* of the version on 7 Sep 2020 [https://github.com/corinalogan/grackles/blob/master/Files/Preregistrations/gspaceuse.Rmd](https://github.com/corinalogan/grackles/blob/master/Files/Preregistrations/gspaceuse.Rmd).

<img width="5%" src="logoOpenAccess.png"> <img width="5%" src="logoOpenCode.png"> <img width="5%" src="logoOpenPeerReview.png">

**This preregistration has been pre-study peer reviewed and received an In Principle Recommendation by:**

Blandine Doligez (2020) Explore and move: a key to success in a changing world? *Peer Community in Ecology*, 100058. [10.24072/pci.ecology.100058](https://doi.org/10.24072/pci.ecology.100058)

 - Reviewers: Joe Nocera, Marion Nicolaus, and Laure Cauchard

### ABSTRACT
Great-tailed grackles (*Quiscalus mexicanus*) are rapidly expanding their geographic range [@wehtje2003range]. Range expansion could be facilitated by consistent behavioural differences between individuals on the range edge and those in other parts of the range [@duckworth2007coupling; @lindstrom2013rapid]. Movement behaviors in particular are thought to be important for expanding and invasive species [@bubb2006movement]. There is evidence for a relationship between individual exploratory traits and dispersal [the permanent movement an individual makes from its birth site to the place where it reproduces; @greenwood1982natal], but it is still unknown whether individual differences in exploration relate to adult daily movement patterns within the home range (“space use”). Evidence suggests that daily space use behavior can vary consistently among individuals [@hertel2020guide], but no studies have examined the relationship between space use behavior and measures of exploration, or space use behavior in different populations along an expanding range. Here we will study grackles from 3 populations that span the extent of the current range - Central America (their original range), Arizona (middle of the northern expanding edge), and northern California (near the northern edge of their range). We will test whether performance on an exploration task in captivity relates to subsequent space use behavior in the wild measured using home range size, the autocorrelation of step length (distance between two sequential observations) and turning angle for each individual over time [@pacheco2019nahua], and the repeatability of each individual's occurrence in particular geographic locations. Results will inform whether individual differences in space use behavior are associated with consistent individual differences in exploration, and whether daily space use patterns differ in populations at different points on an expanding range. If space use behavior correlates with measures of exploration, then space use data could be used to inform conservation management strategies (e.g. identify which individuals are likely to remain in new or restored habitat after a translocation [@may2016predicting]) in species where it is not logistically feasible to *a priori* measure exploration in captivity. If daily space use patterns differ in populations at different locations on the expanding range, then future studies could investigate whether space use and exploration are coupled with other traits such as long-distance natal dispersal to form an “invasion syndrome” [@chapple2012can].

### INTRODUCTION

Problematic range expansions of invasive species are occurring across the globe. However, not all species that invade novel areas outside of their traditional range become established [@chapple2012can; @hayes2008there] and it is still unclear what characteristics facilitate successful invasions [@fogarty2011social]. These characteristics are important for predicting potential invasions, especially since invasive species are implicated as a leading cause of biodiversity loss [@butchart2010global; @clavero2005invasive]. One hypothesis is that the composition of behavioral types might facilitate geographic range expansion [@carere2013animal] because certain individuals may possess traits such as exploration and aggressiveness that make them more likely to succeed when venturing into new habitats and competing with heterospecifics [@cote2010personality]. Evidence is already accumulating that links consistent individual differences in behavioral types and invasion success. For example, @duckworth2007coupling found more aggressive behavior in western bluebird males on the range edge, suggesting that range expansion was facilitated by aggressive males dispersing further and displacing less aggressive mountain bluebirds. Additionally, invasive cane toads on the invasion front in Australia spent longer in the dispersal mode and moved greater distances than toads tracked several years later in the same location [@lindstrom2013rapid]. 

Within-species variation in the ability (movement) and motivation (exploratory tendency) to encounter conspecifics, novel foods, and novel food sources could be a limiting factor in successful species range expansions [@spiegel2016feedback]. In novel areas, the occurrence of conspecifics, food, predators and other environmental factors may not be as easily detectable or recognizable, and may be distributed differently across the landscape than in core areas of the range. Although individuals with exploratory phenotypes may be more successful at colonizing new habitat, exploratory individuals are also at higher risk of predation [@stuber2013slow], and could be less likely to find local food sources [@van2009personality]. Consequently, for establishment in new areas, individuals that exhibit a range of exploratory behavior are needed, and the interaction between space use and exploratory tendency is likely important for finding novel foods and food sources. Additionally, while dispersal [the permanent movement an individual makes from its birth site to the place where it reproduces; @greenwood1982natal] is necessary to initially invade the novel habitat, subsequent daily space use could determine establishment success. For example, on the range edge conspecific density might be lower and individuals may need to use space differently [@bubb2006movement]. However, current research on invasion success and movement behavior focuses heavily on dispersal, and while dispersal and exploratory tendencies have been shown to be associated [@cote2010personality], we do not know how exploratory tendency influences space use patterns in the daily lives of invading individuals.

Space use behavior is expected to be influenced by internal states like exploratory tendency and hunger, as well as the non-random distribution of available habitat and resources [@nathan2008movement]. Space use can also consistently differ among individuals [@hertel2020guide], which indicates that each individual has distinct preferences for how, when, and where to move within its home range. Traditional analyses of animal space use required spatial and temporal independence of data points for statistical analysis [@swihart1985testing], yet movement data are unlikely to meet these criteria. Spatial and temporal autocorrelation occurs when the position of an individual at a given time is tightly linked to its position both before and after, including cases where individuals are repeatedly found in the same locations across time. This autocorrelation is an intrinsic component of space use behavior and eliminating it can reduce biological relevance and obscure relationships with behavioral types [@dray2010exploratory]. Therefore, the autocorrelated nature of movement paths could be important to illuminate the relationship between individual differences in exploratory tendency and daily space use. 

Great-tailed grackles (*Quiscalus mexicanus*, hereafter “grackles”) are rapidly expanding their geographic range [@wehtje2003range]. They are invasive because they meet the criteria for the establishment and spread stages of invasion [@blackburn2011proposed], and they are considered a pest in some areas where they reduce fruit crop yields [@glahn1997controlling]. The nature and level of ecological and social factors grackles experience may vary in importance between populations. Generally, this species is strongly associated with human-modified landscapes and is able to take advantage of a variety of human foods (e.g. crops, at our outdoor cafes, and out of our garbage cans) in addition to foraging on insects and on the ground for natural food items [@johnson2001great]. Furthermore, they exhibit variation in territorial social behavior in both 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 where multiple females place their nests within that territory to raise the young [@johnson2000male]. Roaming males are also present and can obtain extra-pair copulations with females on other males’ territories [@johnson2000male]. 

In this investigation, we aim to understand whether exploratory tendency, measured in captivity, is associated with space use behavior measured in the wild in grackles from three populations that span the current range: Central America (their original range), Arizona (middle of the northern expanding edge), and northern California (near the northern edge of their range). Exploration, measured here following the protocol described in [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]. Previous studies on birds have found a relationship between some movement behaviors (home range size and natal dispersal) in the wild and exploration measured in captivity using novel environment [@minderman2010novel; @dingemanse2003natal] and novel object [@mettke2002significance] tests. Therefore, these measures have the potential to be relevant to grackles. 

After measuring exploratory tendency in captivity, we will release grackles and use telemetry to measure their space use behavior. We will quantify a typical measure of space use that controls for autocorrelation (i.e. home range size), but we will also investigate the behavioral relevance of autocorrelation in space use with two relatively new methods. The first method will describe individual differences in path-level movement behavior by analyzing the autocorrelation of step length (distance between two sequential observations) and turning angle for each individual over time [@pacheco2019nahua; @hertel2020guide], while the second method will describe individual differences in spatial preferences by analyzing the repeatability of each individual's occurrence in particular geographic locations. 

With these data we will test two separate hypotheses. First, we will test whether grackles’ performance on exploration tasks in captivity is related to the space use metrics of the same individuals that we quantify from their movement behavior in the wild. Second, we will test the hypothesis that the space use behavior of wild grackles systematically varies with population location along the extent of the range. These results will have implications for invasive species management [@carere2013animal], but also could be used to inform the conservation of threatened and endangered species [e.g. identify which individuals are likely to remain in new or restored habitat after a translocation; @may2016predicting]. This will be especially relevant for species where it is not logistically feasible to *a priori* measure exploration in captivity. 


### 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). Originally we attached radio tags to grackles released from the aviaries to ensure that we could find their nest sites and track measures of foraging behavior and reproductive success for these individuals for which we have an extensive amount of behavioral data from aviary tests. Now we plan to additionally use the radio tags to collect data for this space use preregistration. This preregistration was written in June 2019, while at the same time increasing the number of GPS points taken per tracking session per bird to provide enough data for the analyses here, and submitted in September 2019 to PCI Ecology for pre-study peer review. Reviews were received in December 2019 and we revised and resubmitted in March 2020. A second round of reviews was received in July 2020 and we revised and resubmitted in August 2020. A third round of reviews was received on September 1, 2020 and we revised and resubmitted on September 15, 2020. A fourth round of reviews was received on September 17, 2020 and we revised and resubmitted on September 30, 2020.

### 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) for methods) 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. Previous studies on birds have found a relationship between movement behavior in the wild and exploration measured in captivity using novel environment [@minderman2010novel; @dingemanse2003natal] and novel object [@mettke2002significance] tests, therefore the measures we are investigating have the potential to be relevant to grackles. We expect space use to vary within an individual across breeding seasons because 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:** More exploratory grackles, i.e. individuals that get closer to the novel object and novel environment will be found in a larger expanse (larger home range size), use less predictable movement patterns (low autocorrelation of step lengths and turning angles), and occupy a greater variety of spatial locations. This would suggest that exploratory individuals are more willing to move into novel areas in the wild.

**Prediction 1 alternative 1:** The more exploratory grackles will be found in a smaller expanse (smaller home range size), use more predictable movement patterns (high autocorrelation of step lengths and turning angles), and consistently occupy the same spatial locations. This would suggest that the more exploratory individuals may dedicate more time to investigating a smaller area within their home range rather than moving 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 would suggest 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 would suggest that, in these populations located in human-modified environments, 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 contained within human-made packaging (e.g. grackles search inside 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 would suggest that the measures of exploration in captivity 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 primarily 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 were novel. 

#### 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). These populations are theoretically connected, however actually moving between two of our field sites within a few grackle lifespans is unlikely due to the large distances between field sites and two geographic barriers (the Sierra Nevada and Sierra Madre mountain ranges, and the high elevation areas of Mexico).

**Prediction 2:** Home range size will increase, autocorrelation of step lengths and turning angles will decrease (i.e. 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, the grackles sampled from our site 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 the sample of grackles in the core of the range (Central America). Grackles in the sample from the middle of the expanding range will be intermediate in space use (Arizona). Such population differences in space use behavior may relate to range expansion because some of the individuals on the leading edge of the range may use more space and move longer distances [@duckworth2007coupling]. However, larger-scale sampling of grackle groups across the strata of the expansion front and core range would be needed to more robustly validate the hypothesis that our cross-site differences are indicative of a broader pattern driven by the location of the expansion front.

**Prediction 2 alternative 1:** Grackles sampled on the edge of the range will have smaller overall home range sizes, high autocorrelation in step length and turning angle (i.e. movement behavior will be more predictable), and consistently use the same spatial locations compared to grackles sampled in the middle or core of the current range. This would suggest that suitable habitat may be distributed in small patches, that novel habitats at the edge of the range may have high predation on grackles that use more space, and/or that individual grackles specialize on certain novel habitat types that are patchily distributed.

**Prediction 2 alternative 2:** We will find no difference across the geographic range in the space use behavior of the grackles sampled. This would suggest that, on average, all grackles may use the same amount of space, or that there is a similar distribution of individual differences in space use in each population. Alternatively, it could indicate that we did not detect differences because we measured adults rather than juveniles. Grackles sampled in different populations may converge on similar space use behavior during development, or juvenile grackles may disperse further on the edge of the range. However, we are not able to detect these differences with our data, which is primarily from adults.

Table 1. Population characteristics for each of the three field sites. The number of generations at a site is based on a generation length of 5.6 years for this species [@GTGRbirdlife2018] and on the first year in which this species was reported to breed at the location (@wehtje2003range for Arizona, Steve Hampton's pers. comm. reported in @pandolfino2009colonization for Woodland, California). The first confirmed nest sighting in Woodland, California was reported in the Yolo Audubon Society's newsletter *The Burrowing Owl* (July 2004), which Steve Hampton shared with Logan. For Central America, there is no data on the first year in which they started breeding because this species originates in this region, therefore we used the age of the species: 800,000 years [@johnsoncicero2004new]. 

```{r table1, eval=TRUE, warning=FALSE, results='asis', echo=TRUE, include=TRUE}
d <- read.csv(url("https://github.com/raw/corinalogan/grackles/master/Files/Preregistrations/gspaceuse_Table1.csv"), header=F, sep=",", stringsAsFactors=F) 
dat<-data.frame(d)

library(reactable)
reactable(dat, highlight=TRUE, bordered=FALSE, compact=TRUE, wrap=TRUE, resizable=TRUE, 
          columns = list(
            V1 = colDef(name="Site"),
            V2 = colDef(name="Range position"),
            V3 = colDef(name="Breeding since"),
            V4 = colDef(name="Number of years breeding"),
            V5 = colDef(name="Average number of generations"),
            V6 = colDef(name="Kilometers to nearest field site"),
            V7 = colDef(name="Nearest field site"),
            V8 = colDef(name="Citation")
          ))
```

### 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 color-marked grackles in this study have one of two different backgrounds: those that do not have radio tags and those that do. First, we opportunistically track color-marked grackles that do not have a radio tag (and thus have not spent time in the aviaries) to compare whether time spent in the aviaries is related to space use behavior. When a color-marked bird is encountered, researchers track it for 20-90 minutes, recording the spatial location every one minute. If the bird goes out of view, researchers attempt to find it again for 15-30 minutes before moving on. Because these data are opportunistic, we do not attempt to balance for sex, but we aim to follow at least 20 non-tagged individuals in each population. 

Second, we applied radio tags to all aviary-tested birds (estimated 20 individuals per population). These subjects are primarily adults who we attempt to balance for sex, and who spent up to six months in an aviary while they participated in behavioral choice tests (see [Logan et al. 2019](http://corinalogan.com/Preregistrations/g_flexmanip.html) for details) and individual differences assays, including measures of exploration in captivity (see [McCune et al. 2019](http://corinalogan.com/Preregistrations/g_exploration.html) for details), as part of other research projects by this lab. For details about the captive environment, please refer to the preregistration associated with this part of the research: [McCune et al. 2019](http://corinalogan.com/Preregistrations/g_exploration.html). Cognitive and behavioral traits are often not fully developed in hatch year birds [i.e. @zucca2007piagetian]. For our aviary cognitive test battery, we avoided the potentially confounding variable of stage of cognitive development by only taking known adult grackles into the aviaries. We identified grackles as adults using eye color, where hatch year birds have brown eyes, but second year and older birds have yellow eyes. While it is possible that cognitive traits continue to develop in second year and older birds, it is impossible to distinguish grackle age after the bird’s first year [@johnson2001great]. 

Before the aviary grackles are released, they are fitted with VHF radio tags (Lotek PipLL (model Ag386), Advanced Telemetry Systems (model A2455) or Holohil Systems Ltd. (model BD-2)) 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, we now use a leg loop harness [methods as in @rappole1991new] made from sutures and secured with crimp beads and cyanoacrylate glue (Vicryl undyed 36in sutures, item number D9389 at eSutures.com; 0.5mm diameter, absorbable so they fall off after one to four months). 

After release, an experimenter will find and follow (each session is called a “track”) each tagged grackle at least four times per week. On each of these tracks, the experimenter follows the focal individual for approximately 1.5 hours, recording a GPS point every one minute, regardless of whether the bird moved [@cushman2005elephants]. Additionally, we aim to balance tracking data equally during morning and afternoon time periods for all grackles. We hoped to track the same individuals across the breeding and nonbreeding seasons. However, the leg loop harnesses often degraded and fell off after 1 - 4 months. Therefore, we have few data points on the same individuals in both seasons and we will instead compare space use in different individuals (rather than within individuals) across seasons. Researchers maintain a distance of at least 30 m and observe the bird with binoculars so the grackle’s behavior is not influenced or artificially changed. If the grackle alarm calls while oriented towards the researcher (indicating the researcher's presence affected the grackle's behavior), all tracking on that individual is stopped for the day. To ensure we capture all locations the individual visits and not just those where they are most easily seen and followed, tagged grackles that move out of sight during tracking are searched for with telemetry until they are found again. 

Exploratory tendency is an intrinsic factor that could be related to space use behaviors.  However, there are also likely alternative variables that may relate to space use behavior in wild grackles that we must control for by including them as covariates in our models. First, we measure energetic condition (described in [Berens et al. 2019](http://corinalogan.com/Preregistrations/gcondition.html)) to account for differences in the physiological mobility that may limit an individual’s space use behaviors [@nathan2008movement]. Secondly, we measure habitat characteristics such as human food sources and available breeding habitat (described in [Logan et al. 2019](http://corinalogan.com/Preregistrations/g_flexforaging.html)) because these factors of the external environment will affect where grackles choose to move or spend time [@nathan2008movement].

Conspecific density has also been shown to affect home range size in other bird species [@flockhart2016factors; @garabedian2018relative]. To control for the possibility that home range size may vary among our populations due to conspecific density rather than exploratory traits, we will use point count surveys to measure grackle population density. We will place 225 point count stations across the landscape encompassing each population (Tempe, AZ; Woodland, CA; Gamboa, Panama). For each study population, the first central point will be randomly placed within 500 m of the center of the study area. The remaining points are placed in a 500 m grid pattern extending out from this central point. In total, the sample area will cover an area that is 7 km by 7 km. Each point will be visited once during the non-breeding season (Sep-Mar). During the survey, researchers will record all grackles visually and aurally detected for six minutes. 
	
#### Sample size rationale
	
We test as many birds as we can during the approximate 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 60 across the three sites (approximately 20 birds per site with the aim that half of the grackles tested at each site are females). We catch grackles with a variety of methods (mist nets, walk-in traps, and bow nets), 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 will primarily track the space use behavior of these ~60 grackles that have radio tags. As described above, we will also opportunistically collect GPS point locations non-tagged grackles to determine whether grackles that were previously in the aviary have different space use behavior from non-aviary-held grackles after their release. We will attempt to match the sample size of aviary birds, and in our Arizona population we currently have over 20 points [the minimum number for reliably calculating home range size; @noonan2019comprehensive] for 33 individuals that have never had radio tags. We aim to acquire more than 20 points on at least 20 non-tagged grackles in the other two populations as well. Additionally, we attach radio tags to birds that are released early because of their lack of willingness to participate in aviary tests (currently 5 individuals) to determine whether space use behavior differs between participatory and non-participatory grackles.
	
#### Data collection stopping rule
	
We will stop collecting GPS location data on tagged and non-tagged birds when home ranges are fully revealed for data collected in both breeding and non-breeding seasons. 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

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.
  - [Point count protocol](https://docs.google.com/document/d/1zDuoI0v7Rv0iwrTMAZSERH7AB4v85qXzxsx5T0daNNo/edit?usp=sharing) for measuring grackle population density in the study area.

#### 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 (see the [separate preregistration](http://corinalogan.com/Preregistrations/g_exploration.html) for details). The time of day that we collect GPS point locations is counterbalanced within and across birds to account for potential 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)

We adapted commonly used methods to test exploratory tendency of the grackles that are temporarily held in aviaries in response to a novel environment and a novel object. Exploratory tendency is measured as the latency to approach and the closest approach distance to a novel object and a novel environment. Exploration assays occur twice for each bird: once near the beginning of their aviary time (“time 1”) and once again approximately 6 weeks later (“time 2”). Habituation may occur between time 1 and time 2, decreasing the novelty of the experimental setup. However, it is common practice to use the same setup across the repeated assays because it is very difficult to predict how threatening a novel object will be to a grackle. Therefore, if we accidentally introduce objects that are much more or much less threatening across the two time periods, this could obscure our ability to determine whether there are consistent individual differences with regard to these particular novel objects. We will analyze whether behavioral responses during assays are repeatable within individuals and whether exploration of a novel environment correlates with exploration of a novel object, indicating they are measures of the same inherent trait. If the two exploration measures are consistent within individuals and correlate with each other, we will choose as the exploration score the variable with the most data. If the two measures do not correlate, we will include both as independent variables.

#### 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 autocorrelation due to the small time period between each of our GPS locations [@noonan2019comprehensive]. This estimate consists of the area enclosing the 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 5 m grid array across the landscape

One model will be run for each dependent variable

#### Independent variables

**P1 and P1 alternatives 1-4**

1) Exploration of novel environment: Latency to approach up 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 up 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)

5) The number of known breeding sites (shade trees, date palms, marsh vegetation [@johnson2001great)] within the home range of each individual (data collected as part of [Logan et al. 2019](http://corinalogan.com/Preregistrations/g_flexforaging.html))

6) The number of human food source areas (dumpsters, cat food bowls, outdoor restaurant seating areas and parking lots) within the home range of each individual (data collected as part of [Logan et al. 2019](http://corinalogan.com/Preregistrations/g_flexforaging.html))

7) Scaled mass index [@peig2009new] as a measure of energetic condition

8) Maximum group size observed across each individual’s focal follows (data collected as part of [Logan et al. 2019](http://corinalogan.com/Preregistrations/g_flexforaging.html))

9) Season (breeding or non-breeding)


**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 expanding edge (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)

4) Population density (number of grackles per square meter in each study area: Arizona, California, Central America)

5)  Season (breeding or non-breeding)

### 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 analyses requiring data from times 1 and 2. Analyses will be conducted in R [current version `r getRversion()`; @rcoreteam] and Stan [version 2.18; @carpenter2017stan]. 

When calculating our dependent variables we are interested in all movements by grackles, therefore we will not exclude any outlier relocations collected during “normal daily activities” [@calenge2011home]. "Normal daily 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 space use movements. 

From the GPS point locations collected on each individual, we will calculate home range size using the akde function in the R packages ctmm [@calabrese2016ctmm] and sf [@pebesma2018simple]. 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 an individual that frequents the same spatial locations across days.

We will quantify region-specific grackle density by fitting a hierarchical model that accounts for imperfect detection. Specifically, we will use the model developed by @amundson2014hierarchical which integrates data on time of detection and distance estimates to account for the probability a bird is available to be detected (pa), and the probability it is detected given it is available (pd), respectively [@nichols2009inferences]. All parameters (density, pa, and pd) will be modeled as a function of region.  Because we expect that vocalization rates will be greater for males than females, we will also model pa as a function of sex. We will extract the estimate of expected point-level density for each region and use this estimate as the covariate in our movement models.

We will then model the relationship between independent variables describing bird-specific data on performance in the exploration tasks and other covariates (as outlined above), and bird-specific movement parameters (e.g. home range size, autocorrelation in step-size and turning angle, repeatability of spatial location preferences) as our dependent variables. We will use linear models and we will ensure assumptions of normality are met by checking that the residuals from our fitted models are normally distributed.

We will additionally test whether time spent in captivity might alter the social behavior of grackles when subsequently released into the wild by testing maximum group size observed across each individual’s focal follows as a function of the individual’s captivity history (the number of days the individual was temporarily held in the aviaries before data collection on space use began).


#### Ability to detect actual effects

To assess the effects that we will be able to detect given the expected sample size, 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). These power analyses are not fully aligned with our study design, and the expected effect sizes are difficult to estimate due to the lack of prior data on this species; yet we are unaware of better options at this time.

#### 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 20 or more relocations
### Count the number of relocations for each bird, then exclude individuals with less than 20
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<20),]
pts.min$Bird.Name = as.factor(as.character(pts.min$Bird.Name))

#### Calculating home range size ####
#Turn dataframe into telemetry object
telem = pts.min[,c(1,2,4,5,7)] #only need bird name, date, lat and long, time
telem$timestamp = paste(telem$Date,telem$Time, sep=' ') #telem format requires date/time merged together for timestamp
telem = telem[,-c(2,5)] #remove other date and time columns
colnames(telem) = c("individual.local.identifier","location.lat","location.long","timestamp") #use required column names
telem = as.telemetry(telem, timeformat = "%Y-%m-%d %H:%M:%S",timezone = "")

for(i in 1:length(telem)){
  GUESS = ctmm.guess(telem$i, interactive = FALSE)
  m = ctmm.fit(telem$i, GUESS)
  hr = akde(telem$i, m)
  est[i] = summary(hr)$CI[2]
  results = rbind(results,est[i])
} 


### Normalize data for repeatability and regression models 
hist(log(hr$area)) 

### Repeatability of home range in breeding and non-breeding seasons
hr_rpt = rpt(log(area) ~ (1|BirdID), grname = "BirdID", data = hr, datatype = "Gaussian", nboot = 500, npermut = 500)
summary(hr_rpt)

### 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=60). The protocol of the power analysis is here:

*Input:*

Effect size f²                	=	0,24

α err prob                    	=	0,05

Power (1-β err prob)          	=	0,7

Number of predictors          	=	8

*Output:*

Noncentrality parameter λ     	=	14,4000000

Critical F                    	=	2,1260234

Numerator df                  	=	8

Denominator df                	=	51

Total sample size             	=	60

Actual power                  	=	0,7039242

This means that, with our minimum sample size of 60, 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 + History + Breeding + Feeding + Condition + Group + Season, data = data)
hist(m1$resid)
summary(m1)

# Step length 
#m2 = lm(log(std_step)) ~ ExpObj + ExpEnv + Sex + History + Breeding + Feeding + Condition + Group + Season, data = data)
hist(m2$resid)
summary(m2)

# Turning angle
#m3 = lm(log(std_angle)) ~ ExpObj + ExpEnv + Sex + History + Breeding + Feeding + Condition + Group + Season, data = data)
hist(H3$resid)
summary(m3)

# Spatial preferences
#m4 = lm(log(loc_pref)) ~ ExpObj + ExpEnv + Sex + History + Breeding + Feeding + Condition + Group + Season, 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,18

α err prob                    	=	0,05

Power (1-β err prob)          	=	0,7

Number of predictors          	=	4

*Output:*

Noncentrality parameter λ     	=	10,6200000

Critical F                    	=	2,5429175

Numerator df                  	=	4

Denominator df                	=	54

Total sample size             	=	59

Actual power                  	=	0,7025819

This means that, with our minimum sample size of 60, 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 + History + Season, data = data)
hist(m1$resid)
summary(m1)

# Step length
#m2 = lm(log(std_step)) ~ Site + Sex + History + Season, data = data)
hist(m2$resid)
summary(m2)

# Turning angle
#m3 = lm(log(std_angle)) ~ Site + Sex + History + Season, data = data)
hist(m3$resid)
summary(m3)

# Spatial preference
#m4 = lm(log(loc_pref)) ~ Site + Sex + History + Season, 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], SP639866 [2019], and SP402153 [2020])
4) Institutional Animal Care and Use Committee at Arizona State University (protocol number 17-1594R)

### F. AUTHOR CONTRIBUTIONS

**McCune:** Hypothesis development, data collection (trapping, GPS tracking), data analysis and interpretation, write up, revising/editing.

**Folsom:** Data collection (trapping, GPS tracking), revising/editing.

**Ross:** Model development, data analysis and interpretation, revising/editing.

**Bergeron:** Data collection (trapping, GPS tracking), 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 the local PI on the ASU IACUC; Kristine Johnson for technical advice on great-tailed grackles; Julia Cissewski and Sophie Kaube 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, Amanda Overholt, Emily Blackwell, Kaylee Delcid, Brynna Hood, Sam Bowser, Sierra Planck, Samuel Munoz, Sawyer Lung.

### K. [REFERENCES](MyLibrary.bib)