NetworksSummer2021-IndividualMovement.Rmd

---
title: "Individual Movement Analysis"
author: "Matina Donaldson-Matasci"
date: "10/11/2021"
output:
  pdf_document: 
      latex_engine: xelatex
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
require(ggplot2)
require(dplyr)
require(tidyr)
require(ggthemes)
require(grDevices)
require(lme4)
require(glmmTMB)
require(bbmle) #for AICtab

# theme_set(theme_test(base_size=18)) # For poster
theme_set(theme_test(base_size=10)) # For paper

nest.palette <- palette.colors(n = 9, palette = "Classic Tableau")[c(1:7,9)]
nest.names <- c("A2","B3","C1","D2","E1","F1","G1","H1")
names(nest.palette) = nest.names
```

## Hypotheses / Questions
We would like to discover the rules by which individual turtle ants move through a branching structure.

1. Are they more likely to take left or right turns?
1. Are they more likely to take a wider, less angled turn or a narrower, more angled turn?
1. Are they more likely to move inwards towards the trunk, or outwards towards the tips?

We also want to understand how these choices at individual junctions produce paths to specific nests.

1. Which nests are they most likely to visit?
1. Which nests do they reach most quickly?

## Data description

We created a branching structure for the ants to navigate, consisting of a series of many Y-junctions. Each Y-junction has a thicker, main branch leading back to the trunk, and two alternate choices, one of which has a small angle of deviation (10 degrees) and the same diameter as the main branch, and the other of which has a large angle of deviation (70 degrees) and a smaller diameter.

First, given a table of the junction labels and their connections to one another, we would like to construct a table of the different types of turns an ant can take. 

Here is the table of junction entry/exit points, which we will call branches. We are only interested in the columns "Junction" (the junction label), "Direction" (whether it is the main, left or right branch), "Width" (the diameter of the branch) and "Leads.to.Node" (which node the branch leads to). We check to see that it has the right shape: we expect 31 junctions, each of which has a main, left and right branch. So there are 93 rows, one row for each branch of each junction. There are 4 different branch diameters. Most junctions have two different branch widths, but there is one junction that has only 1. 

```{r}
branches <- read.csv("data/NetworksSummer2021-JunctionReference.csv",
                     colClasses=c("Junction"="character",
                                  "Direction"="factor",
                                  "Width"="numeric",
                                  "Leads.to.Node"="character"),
                     na.strings="")
branches <- branches %>%
  dplyr::select(Junction, Direction, Width, Leads.to.Node)
head(branches)
branches %>% summarize(num_junctions=n_distinct(Junction))
branches %>% group_by(Junction) %>% summarize(num_branches=n())
branches %>% group_by(Width) %>% summarize(num_branches=n())
branches %>% group_by(Junction) %>% summarize(num_widths=n_distinct(Width))
```

Now we want to create a table of turns. This consists of a sequence of two things: entering a junction from one branch, and exiting it along (potentially) another branch. That is, for each junction, we want all combinations of branches. To do this, we use a full join, based on the junction. For each turn, we then specify whether it is a left turn or a right turn, based on the entry and exit direction.

This action is then carried out with narrow and wide junction decisions as well.

```{r}
branches <- branches %>%
  group_by(Junction) %>%
  mutate(Rel.Width = if_else(Width > min(Width),"wide","narrow")) %>%
  ungroup()

turns <- full_join(branches, branches,
                   by=c("Junction"), suffix=c(".From",".To"),
                   relationship = "many-to-many") %>%
  rename(Node.From=Leads.to.Node.From, Node.To=Leads.to.Node.To)
 
turns <- turns %>%
  mutate(Turn.Type = factor(case_when(
    Direction.From=="main" & Direction.To=="left" ~ "L",
    Direction.From=="main" & Direction.To=="right" ~ "R",
    Direction.From=="left" & Direction.To=="right" ~ "L",
    Direction.From=="left" & Direction.To=="main" ~ "R",
    Direction.From=="right" & Direction.To=="main" ~ "L",
    Direction.From=="right" & Direction.To=="left" ~ "R",
    TRUE ~ "U"
  )))

turns <- turns %>%
  mutate(Turn.Width = if_else(Node.From == Node.To, NA_character_,case_when(
    Rel.Width.From=="wide" & Rel.Width.To=="narrow" ~ "N",
    Rel.Width.From=="wide" & Rel.Width.To=="wide" ~ "W",
    Rel.Width.From=="narrow" ~ "=")))

turns <- turns %>%
  mutate(Turn.Trunk = if_else(Direction.From == Direction.To, NA_character_,case_when(
    Direction.From=="left" & Direction.To=="main" ~ "to.trunk",
    Direction.From=="right" & Direction.To=="main" ~ "to.trunk",
    Direction.From=="main" & Direction.To=="left" ~ "away.from.trunk",
    Direction.From=="main" & Direction.To=="right" ~ "away.from.trunk",
    Direction.From=="left" & Direction.To=="right" ~ "away.from.trunk",
    Direction.From=="right" & Direction.To=="left" ~ "away.from.trunk"
  )))

```

Finally, we want to join our data on individual movement between junctions with this data describing the types of turns. To do this, first we need to fill in some missing information: where the ant came from, and where it's going. We begin with information about the sequence of junctions visited ("Junction"). This is typically sufficient to reconstruct the whole trajectory: the ant came from the previous junction in the list, and went towards to next junction in the list. However, if the ant made a U-turn along a branch, we have an exception: the ant actually never reached the node it was going towards when it first exited the junction, which is also the node it was traveling from when it reenters the same junction again. In this case, we recorded the junction it was traveling towards in the column ("Exited Toward").

```{r}
movement <- read.csv("data/NetworksSummer2021-IndividualMovement.csv",
                     colClasses="character", na.strings="")
movement <- movement %>%
  dplyr::select(Date, Colony, Ant.ID, Time, Junction, 
         Entered.From, Exited.Toward, Action..inspect.ignore.) %>%
  mutate(Colony=factor(Colony),Ant.ID=factor(Ant.ID)) %>%
  group_by(Colony,Ant.ID) %>%
  # Ants come in pairs, labeled e.g. as 1A and 1B
  # We want to be able to distinguish the first ant (A) 
  # from the second (B) in each pair
  separate(Ant.ID, into=c("Pair.ID","Ant.In.Pair"), 
           sep=1, remove=F, convert=F) %>%
  mutate(Pair.ID=as.factor(Pair.ID),
         Ant.In.Pair=as.factor(Ant.In.Pair)) %>%
  # Here we count up how many times an ant has inspected food
  # or a nest so we can distinguish outgoing & returning paths
  mutate(Inspect=if_else(grepl("inspect",
                               Action..inspect.ignore.),T,F)) %>%
  mutate(Num.Inspections=cumsum(Inspect)) %>%
  # Here we infer the direction an ant came from by 
  # the prior junction it visited and where it went by the
  # next junction it visited, unless that information was given
  # in the Exited.Toward column
  mutate(Node.From=lag(Junction),Node.To=lead(Junction)) %>%
  mutate(Node.From=replace_na(Node.From,"-1")) %>% 
  mutate(Node.To=if_else(is.na(Exited.Toward),Node.To,Exited.Toward)) %>%
  mutate(Node.To=replace_na(Node.To,"-1")) %>% 
  mutate(Entered.From=lag(Exited.Toward))  %>%
  mutate(Node.From=if_else(is.na(Entered.From),Node.From,Entered.From))
```

Now, finally, we should be able to join the two tables together to generate a list of all the turns made by each ant.

```{r}
turn_choices <- left_join(movement,turns,by=c("Junction","Node.From","Node.To")) %>%
  filter(Junction %in% unique(branches$Junction))
write.csv(turn_choices,file="data/NetworksSummer2021-TurnChoices.csv")
```

Besides the turns they made at each junction, we also want to analyze which tips (ends of branches, containing food or nests) they reached. To figure this out, we look for visits to junctions with labels that are in our list of tips. For each such visit, we ask how many nodes they passed along the way since the start (a measure of path length).

```{r}
tips <- read.csv("data/NetworksSummer2021-TipReference.csv",
                 colClasses=c("Tip"="character","Type"="factor",
                              "Location"="factor",
                              "Distance.from.Trunk"="numeric"),
                 na.strings="")

movement <- movement %>%
  group_by(Colony,Ant.ID) %>%
  mutate(Num.Nodes.Passed=row_number())

tip_visits <- inner_join(movement,tips,by=c("Junction"="Tip")) %>%
  rename("Tip.Location"="Location")

write.csv(tip_visits,file="data/NetworksSummer2021-TipVisits.csv")
```

## Data analysis

### Which nests are they most likely to visit?

Of the eight different nests, which one are the ants most likely to visit first? We filter out only visits to nests, and then take the first one for each ant. The plot shows that many ants visited nests A2 and B3 first.

```{r}
first_nest_visits <- tip_visits %>%
  filter(Type=="nest") %>%
  group_by(Colony,Ant.ID) %>%
  summarize(Nest.Location=factor(first(Tip.Location),levels=nest.names),
            Junction=first(Junction),
            Num.Nodes.Passed=first(Num.Nodes.Passed))

ggplot(first_nest_visits, aes(x=Nest.Location,fill=Nest.Location)) + 
  geom_bar() + ylim(0,14) +
  scale_fill_manual(values=nest.palette,drop=F) +
  scale_x_discrete(drop=F) +
  xlab("Nest Location") + 
  ylab("Ants reaching this nest first") +
#  labs(title = "Nest Discovery (Individual Experiments)") +
  theme(legend.position="none",plot.title = element_text(hjust = 0.5))
#ggsave("images/FirstNestVisits-NumAnts-POSTER.png", 
#       dpi=300, width = 8, height = 5, units = "in")
ggsave("images/FirstNestVisits-NumAnts-PAPER.png", 
       dpi=600, width = 3, height = 2.5, units = "in")
```

Is there a significant preference for particular nests over others? We can do a chi-squared goodness of fit test to test the null hypothesis that ants are equally likely to reach all 8 nests first. Assumptions are: (1) independence of observations, and (2) no more than 20% of the expected values in each cell are less than 5. Here we have a total of 46 observations in 8 categories, so the expected value for each cell is 5.75.

```{r}
chisq.test(table(first_nest_visits$Nest.Location))
```

This shows that there is a significant preference for certain nests over others (chi-squared goodness of fit, X2=27.7, df=7, p=0.0002).

We also would like to know how direct their route to get there was. If they take a direct route, this will be easier to reinforce with pheromone. To figure this out, we plot the number of nodes each ant passed along the way to the first nest they visited. The shortest possible path is 6 nodes long, because that is the distance from the entry point to each tip.

```{r}
ggplot(first_nest_visits, aes(x=Nest.Location,
                              y=Num.Nodes.Passed)) + 
  geom_boxplot(aes(fill=Nest.Location), outlier.shape=NA) + 
  geom_jitter(height=0,width=0.2, shape=21) +
  scale_fill_manual(values=nest.palette,drop=F) +
  scale_x_discrete(drop=F) +
  geom_abline(slope=0,intercept=6,lty=2) +
  xlab("Nest Number") + ylab("Path length (number of nodes)")  +
#  labs(title = "First Nest Path Length (Individual Experiments)") +
  theme(legend.position="none", plot.title = element_text(hjust = 0.5))
# ggsave("images/FirstNestVisits-PathLength.png", 
#        dpi=300, width = 8, height = 5, units = "in")
ggsave("images/FirstNestVisits-PathLength-PAPER.png", 
       dpi=600, width = 3, height = 2.5, units = "in")
```

The plot suggests that nest B3 has the most direct routes; in fact there are several ants that went directly there -- something that never happened for any other nest. 

```{r}
minNodes <- first_nest_visits %>%
  ungroup() %>%
  summarize(minNodes=min(Num.Nodes.Passed)) %>%
  pull(minNodes)

first_nest_visits %>% 
  group_by(Nest.Location) %>%
  filter(Num.Nodes.Passed == minNodes)
```

Let's test this using a negative binomial model to represent the number of nodes an ant passed through on the way to its first nest. We'll use 'sum contrasts' for the nest location, which means we'll compare the first path length for each nest to the mean first path length observed across all nests.

```{r}
model.pathlength.bynest <- glmmTMB(Num.Nodes.Passed ~ Nest.Location + 
                          (1|Colony),
                          contrasts = list(Nest.Location = "contr.sum"),
                        family="nbinom2", data=first_nest_visits)

sjPlot::tab_model(model.pathlength.bynest)
```

This shows that the overall path length is 22.52 nodes; the first nest location (A2) has significantly longer paths (39% more nodes), the second nest location (B3) has significantly shorter paths (30% fewer nodes), and the third nest location (D2) has significantly longer paths (95% more nodes). These three are likely the only ones that are significant because the other nests received fewer visits (5 or less).

### Are individual ants' nest choices influenced by the previous ant in the pair?

The previous analysis assumes that all ants' choices are independent, but the two ants in a pair may not be. Here we determine whether the second ant in a pair is more likely to choose the same nest as the first ant in a pair than one would expect by chance.

```{r}
paired_first_nest_visits <- first_nest_visits %>% 
  separate_wider_position(Ant.ID,
                    widths=c(Pair.ID=1,
                             Pair.Order=1)) %>%
  pivot_wider(id_cols=c(Colony,Pair.ID), 
              names_from = Pair.Order, 
              values_from= 
                c(Nest.Location,Num.Nodes.Passed))
```

### Do ants have a turning angle bias?

We would like to know whether ants choose which direction to turn based on the turning angle. In particular, we hypothesize that they prefer a shallower (less deviating) turn at an asymmetric Y-shaped junction. 

To determine this, we will ask whether ants tend to take the sharper (more deviating) or shallower (less deviating) turn, and whether any observed preference depends on:

1. Whether the sharper turn is on the left side or the right side

1. The direction the ant approaches the junction from: the main branch (main), the thinner more angled offshoot (sharp), or the thicker less angled offshoot (shallow)

To do this, we'll need to introduce some new variables to describe the junctions and the turns at each junction. First of all, we classify each junction as a left-handed one (with the narrower, sharper-angled offshoot to the left) or a right-handed one (with the narrower, sharper-angled offshoot to the right). There are two special cases where the width is the same for all three branches at the junction, but the difference in angle nonetheless gives them a handedness.

```{r}
junction_type <- branches %>%
  filter(Rel.Width=="narrow") %>%
  mutate(Handedness_raw = case_when(
    Direction=="left" ~ "LH", #left-handed junction
    Direction=="right" ~ "RH", #right-handed junction
    Junction == 111 ~ "RH", #can't tell from width, but it's RH
    Junction == 221 ~ "LH" # can't tell from width, but it's LH
  )) %>%
  group_by(Junction) %>%
  summarize(Junction.Handedness = 
              names(which.max(table(Handedness_raw)))) %>%
  mutate(Junction.Handedness=factor(Junction.Handedness))
```

Now we can create a new way to classify each turn (`Turn.Angle`), as "sharp" or "shallow", according to how the ant entered the junction (`Direction.From`: "main", "left", "right") and the handedness of that junction (`Junction.Handedness`: "LH", "RH"). For example, if the ant enters from "main", and turns left on a left-handed junction, that means it has taken the sharper turn.

We also classify the direction it has come from (`Angle.From`) according to the angle, so for example an ant coming from the left branch on a left-handed junction is coming from the "sharp" direction (the more-angled offshoot), while an ant coming from the left branch on a right-handed junction is coming from the "shallow" direction (the less-angled offshoot).

Finally, we determine whether the sharper turn is on the left side or the right side (`sharpIsLeft`), which depends on whether it is a left or right-handed junction, and which branch the ant is coming from.

```{r}
turns <- turns %>%
  left_join(junction_type, by="Junction") 

turns <- turns %>%
  mutate(
    Turn.Angle = case_when(
      Direction.From=='main' & Turn.Type=='L' & Junction.Handedness=="LH" ~ "sharp",
      Direction.From=='main' & Turn.Type=='R' & Junction.Handedness=="LH" ~ "shallow",
      Direction.From=='main' & Turn.Type=='L' & Junction.Handedness=="RH" ~ "shallow",
      Direction.From=='main' & Turn.Type=='R' & Junction.Handedness=="RH" ~ "sharp",
      Direction.From=='left' & Turn.Type=='L' ~ "sharp",
      Direction.From=='left' & Turn.Type=='R' ~ "shallow",
      Direction.From=='right' & Turn.Type=='L' ~ "shallow",
      Direction.From=='right' & Turn.Type=='R' ~ "sharp"),
      #Junction==111 & Turn.Type=='L' ~ "shallow",
      #Junction==111 & Turn.Type=='R' ~ "sharp",
      #Junction==221 & Turn.Type=='L' ~ "sharp",
      #Junction==221 & Turn.Type=='R' ~ "shallow")) %>%
    isUTurn = Turn.Type=="U",
    isSharpTurn = Turn.Angle == "sharp",
    sharpIsLeft = case_when(
      Direction.From=='main' & Junction.Handedness=="LH" ~ TRUE,
      Direction.From=='left' & Junction.Handedness=="LH" ~ TRUE,
      Direction.From=='left' & Junction.Handedness=="RH" ~ TRUE,
      TRUE ~ FALSE # for all other cases, sharp is not left
    ),
    Angle.From = factor(case_when(
      Direction.From == 'main' ~ "main",
      Direction.From == 'left' & Junction.Handedness == "LH" ~ "sharp",
      Direction.From == 'left' & Junction.Handedness == "RH" ~ "shallow",
      Direction.From == 'right' & Junction.Handedness == "LH" ~ "shallow",
      Direction.From == 'right' & Junction.Handedness == "RH" ~ "sharp"
      # Direction.From == 'left' & Junction == 111 ~ "shallow",
      # Direction.From == 'right' & Junction == 111 ~ "sharp",
      # Direction.From == 'left' & Junction == 221 ~ "sharp",
      # Direction.From == 'right' & Junction == 221 ~ "shallow"
    ))
  )
```

Finally we join the actual turn choices individual ants made (`turn_choices`) with the different types of turns (`turns`), using `Node.From` and `Node.To`. This creates a table of turns taken, including our new columns (`Turn.Angle`, `Junction.Handedness`, `Angle.From`, `sharpIsLeft` and `isSharpTurn`). For this analysis we only use exploratory travel (all turns before the first nest is reached) by the first ant in the pair, and we ignore U-turns.

```{r}
turn_choices_angle <- turn_choices %>%
  left_join(dplyr::select(
    turns, Node.From, Node.To, Turn.Angle, Angle.From,
    isSharpTurn, sharpIsLeft, Junction.Handedness), 
    by=c("Node.From", "Node.To")) %>%
  filter(Num.Inspections==0,Ant.In.Pair=="A") %>%
  filter(Turn.Type != "U")
```

Now we can analyze the turn choices made. We use logistic regression (a binomial generalized linear mixed effects model) to create a model of whether the ant chooses the shallower or sharper turn at each junction. To begin with we look at the null model, which tells us overall whether the turn choice depends on angle at all.

```{r}
null_model <- glmer(
  isSharpTurn ~ (1|Colony:Ant.ID), 
  family = binomial, turn_choices_angle)
summary(null_model)
```

The fact that the intercept is significantly below zero means that ants are significantly less likely to take the sharper turn than the shallower turn. Now let's see if that bias is influenced by anything else. First we can see whether those sharp turns are more likely to be taken if they are also on the left side.

```{r}
leftright_model <- glmer(
  isSharpTurn ~ sharpIsLeft + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_angle)
summary(leftright_model)
anova(null_model,leftright_model)
```

The answer is yes; if the sharp turn is on the left, it is significantly less likely to be taken. Now let's see whether this depends on the direction it came from (`Angle.From`: the main branch, the sharper-angled offshoot, or the shallower-angled offshoot). We'll construct a series of nested models including the `Angle.From` alone, including `Angle.From` and `sharpIsLeft` as additive effects, and including all combinations of `Angle.From` and `sharpIsLeft`. We can compare models using `anova` (a likelihood ratio test) to see whether a more complex model provides a significantly better fit.

```{r}
# We'll try three successively more complex models
angle_model <- glmer(
  isSharpTurn ~ Angle.From + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_angle)
angle_model_plus <- glmer(
  isSharpTurn ~ Angle.From + sharpIsLeft + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_angle)
angle_model_star <- glmer(
  isSharpTurn ~ Angle.From * sharpIsLeft + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_angle)

# Which model is best?
AICtab(null_model, leftright_model,
       angle_model,angle_model_plus,angle_model_star)
```


```{r}
summary(angle_model_plus)
```

The best model is the one including both the direction of approach, and whether the sharp turn is on the left side, as independent additive effects on the log-odds ratio. When the ant is approaching the turn from the less-angled offshoot, it is much more likely to continue straight back along the main branch than to turn sharply onto the other offshoot. 

To interpret the model, let's look at the estimates and 95% CI for the coefficients for each fixed effect (which direction the ant is coming from, and whether the sharp turn is on the left). We can also use the model to make predictions for each possible approach.

```{r}
# these confidence intervals are in the log-odds scale
#angle.ci <- confint(angle_model_plus,parm="beta_")
#angle.ci

# create a small table containing all six types of approaches
approach.types <- turn_choices_angle %>% 
  ungroup() %>% 
  dplyr::select(Angle.From, Junction.Handedness, Direction.From, sharpIsLeft) %>% 
  distinct() %>%
  arrange(Angle.From, sharpIsLeft)
# use the table to make predictions from the model
approach.types$Sharp.Turn.Prob <- predict(
  angle_model_plus, newdata=approach.types, 
  re.form=NA, # this says ignore random effects in the predictions
  type="response" # this says convert output to probabilities
)
# approach.types <- approach.types %>%
#   mutate(Angle.From = case_when(
#            Angle.From == "main" ~ "Main",
#            Angle.From == "sharp" ~ "Secondary",
#            Angle.From == "shallow" ~ "Primary"))
```

## Visualization of turning angle bias data

Here are a couple of ideas for visualization (not mutually exclusive):

1. Do a plot with all turns (excluding U-turns) on the x axis, and sharp turns on the y axis, broken out into six panels by `Angle.From` and `Junction.Handedness`. Each point represents one junction (node), and you could use colors and/or shapes to distinguish among junctions. The predictions above give you a regression line to plot for each panel.

2. For each of the six combinations of `Angle.From` and `Junction.Handedness`, draw colored arrows (approach, sharp turn, shallow turn) on top of a diagram of the junction, showing the direction of approach and each of the two turns. The width of each arrow should be proportional to the number of ants that approached or turned, summed across all junctions in the category and you can write the actual number next to it or inside it.

Here we'll create the first plot. 


```{r}
angle_graph_data <- turn_choices_angle %>%
  ungroup() %>%
  group_by(Angle.From, Junction.Handedness, Junction) %>%
  dplyr::select(isSharpTurn) %>%
  filter(Junction != 0) %>%
  summarize(Total = n(), Sharp = sum(isSharpTurn)) %>%
  mutate(Junction = as.numeric(Junction),
         Level = case_when(
           Junction < 10 ~ "1",
           Junction >= 10 & Junction < 100 ~ "2",
           Junction >= 100 & Junction < 1000 ~ "3",
           Junction >= 1000 ~ "4"),
         LogTotal = log(Total + 1),
         LogSharp = log(Sharp + 1)#,
         # Angle.From = case_when(
         #   Angle.From == "main" ~ "Main",
         #   Angle.From == "sharp" ~ "Secondary",
         #   Angle.From == "shallow" ~ "Primary"
         # )
    )

fixed_effects <- 
  as_tibble(effects::Effect(
    focal.predictors=c("Angle.From",
        "sharpIsLeft"), 
    mod=angle_model_plus)) #%>%
  # mutate(Angle.From = case_when(
  #          Angle.From == "main" ~ "Main",
  #          Angle.From == "sharp" ~ "Secondary",
  #          Angle.From == "shallow" ~ "Primary"))

angle_graph_data
fixed_effects
```

Join the fixed_effects data frame with approach.types to get them to match for more variable usage.

```{r}
approach.types <- approach.types
fixed_effects <- as_tibble(fixed_effects) %>% 
  mutate(sharpIsLeft = as.logical(sharpIsLeft))
approach.types
fixed_effects

angle_fit_data <- left_join(approach.types, fixed_effects, by = c("Angle.From","sharpIsLeft")) #%>%
   # mutate(Angle.From = case_when(
   #          Angle.From == "main" ~ "Main",
   #          Angle.From == "sharp" ~ "Secondary",
   #          Angle.From == "shallow" ~ "Primary"))
angle_fit_data
```

```{r}
sjPlot:: tab_model(angle_model_plus)
```

Plot the effects of angle_model_plus

```{r}
# Use the effects package --> effect function. term= the fixed effect you want to get data on, mod= name of your model.

#effects_angle <- effects::effect(term= "Angle.From", mod= angle_model_plus)
#summary(effects_angle) #output of what the values are
#df_angle <- as.data.frame(effects_angle)
```

Non-transformed data visualization:

```{r}
maxtotal <- max(angle_graph_data$Total)
maxsharp <- max(angle_graph_data$Sharp)
maxtrue <- max(maxtotal,maxsharp*2)

testpolygon2 <- tibble(xu = c(0,0,maxtrue),
                       yu = c(0,maxtrue/2,maxtrue/2),
                       xl = c(0,maxtrue,maxtrue),
                       yl = c(0,0,maxtrue/2))

angle_ribbon_data_pre <- angle_fit_data %>%
  mutate(yminf = lower * maxtotal, ymaxf = upper * maxtotal, x0 = 0, xf = maxtotal, ymin0 = 0, ymax0 = 0)

angle_ribbon_data_x <- angle_ribbon_data_pre %>%
  pivot_longer(cols = c(x0,xf), values_to = "x") %>%
  dplyr::select(Junction.Handedness, Angle.From, x)

angle_ribbon_data_ymax <- angle_ribbon_data_pre %>%
  pivot_longer(cols = c(ymax0,ymaxf), values_to = "ymax") %>%
  dplyr::select(ymax)

angle_ribbon_data_ymin <- angle_ribbon_data_pre %>%
  pivot_longer(cols = c(ymin0,yminf), values_to = "ymin") %>%
  dplyr::select(ymin)

angle_ribbon_data <- bind_cols(angle_ribbon_data_x,angle_ribbon_data_ymax) %>%
  bind_cols(angle_ribbon_data_ymin)

angle_ribbon_data_pre2 <- angle_fit_data %>%
  mutate(yminf = lower * maxtrue, ymaxf = upper * maxtrue, x0 = 0, xf = maxtrue, ymin0 = 0, ymax0 = 0)

angle_ribbon_data_x2 <- angle_ribbon_data_pre2 %>%
  pivot_longer(cols = c(x0,xf), values_to = "x") %>%
  dplyr::select(Junction.Handedness, Angle.From, x)

angle_ribbon_data_ymax2 <- angle_ribbon_data_pre2 %>%
  pivot_longer(cols = c(ymax0,ymaxf), values_to = "ymax") %>%
  dplyr::select(ymax)

angle_ribbon_data_ymin2 <- angle_ribbon_data_pre2 %>%
  pivot_longer(cols = c(ymin0,yminf), values_to = "ymin") %>%
  dplyr::select(ymin)

angle_ribbon_data2 <- bind_cols(angle_ribbon_data_x2,angle_ribbon_data_ymax2) %>%
  bind_cols(angle_ribbon_data_ymin2)
```


```{r}
junction.labels <- c(LH="Left-handed", RH="Right-handed")
angle.labels <- c(main="From main",shallow="From primary",sharp="From secondary")
angle_graph_data %>%
  ggplot() +
  facet_grid(cols = vars(Junction.Handedness), rows = vars(Angle.From),
             labeller=labeller(Junction.Handedness=junction.labels,
                               Angle.From=angle.labels)) +
  geom_point(mapping = aes(x = Total, y = Sharp)) +
  geom_abline(data=angle_fit_data, aes(slope = fit, intercept = 0),size = 1) +
  geom_ribbon(data=angle_ribbon_data2, aes(x = x, ymin = ymin, ymax = ymax), alpha = .3) +
  geom_polygon(data = testpolygon2, aes(x = xu, y = yu), alpha = .25, fill = "red") +
  geom_polygon(data = testpolygon2, aes(x = xl, y = yl), alpha = .25, fill = "blue") +
  coord_cartesian(xlim = c(0,maxtrue), ylim = c(0, maxtrue/2), 
                  expand = FALSE, clip="on") +
  labs(y = "Number of Sharp Turns", x = "Number of Total Turns")
ggsave("images/Angle_Data_Visualization-PAPER.png",dpi=600, width = 3.25, height = 4.5, units = "in")
```

Here is a table of the total turns in each category to be used in creating the second figure.

```{r}
turn_choices_angle %>%
  ungroup() %>%
  group_by(Direction.From, Junction.Handedness) %>%
  dplyr::select(isSharpTurn) %>%
  summarize(Total = n(), Sharp = sum(isSharpTurn))
```

### Do ants have U-turn bias?

We would like to know whether ants choose whether to make a U-turn based on the turning angles available to them. If ants are more likely to take a U-turn at certain types of junctions, coming from certain directions, this could induce an overall bias in their paths.

To determine this, we will ask how likely it is for ants to make a U-turn overall, and whether this probability depends on:

1. Whether the sharper turn is on the left side or the right side

1. The direction the ant approaches the junction from: the main branch (main), the thinner more angled offshoot (sharp), or the thicker less angled offshoot (shallow)

That is, we will use exactly the same explanatory variables as in the turning angle model, but here change the response variable to be whether they made a U-turn or not. First we create the data frame with all the additional variables about the turning context from before, but this time we include the U-turns.

```{r}
turn_choices_u <- turn_choices %>%
  left_join(dplyr::select(
    turns, Node.From, Node.To, Turn.Angle, Angle.From, isUTurn, sharpIsLeft, Junction.Handedness), 
    by=c("Node.From", "Node.To")) %>%
  filter(Num.Inspections==0,Ant.In.Pair=="A")
```

Now we can analyze the turn choices made. We use logistic regression (a binomial generalized linear mixed effects model) to create a model of whether the ant makes a U-turn at each junction. To begin with we look at the null model, which tells us overall what the chance of making a U-turn is.

```{r}
null_model_u <- glmer(
  isUTurn ~ (1|Colony:Ant.ID), 
  family = binomial, turn_choices_u)
sjPlot::tab_model(null_model_u)
```

This shows that the odds ratio of making a U-turn is about 0.22 across all junction types (out of 122 turns, 22 will be U-turns). Now let's see if that ratio is influenced by anything else. We'll try various combinations of (1) whether the sharp turn is on the left (`sharpIsLeft`), and (2) whether the ant is approaching from the main branch (`Angle.From`="main"), the primary branch, which is thicker and less angled away from the main (`Angle.From`="shallow"), or the secondary branch, which is thinner and more angled away from the main (`Angle.From`="shallow"). We can compare models using AIC (the Akaike Information Criterion) to decide which one best explains the observations without fitting too many unnecessary parameters.

```{r}
# We'll try several models

leftright_model_u <- glmer(
  isUTurn ~ sharpIsLeft + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_u)
angle_model_u <- glmer(
  isUTurn ~ Angle.From + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_u)
angle_model_plus_u <- glmer(
  isUTurn ~ Angle.From + sharpIsLeft + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_u)
angle_model_star_u <- glmer(
  isUTurn ~ Angle.From * sharpIsLeft + (1|Colony:Ant.ID), 
  family = binomial, turn_choices_u)

# Which model is best?
AICtab(null_model_u,
       angle_model_u,
       leftright_model_u,
       angle_model_plus_u,
       angle_model_star_u)

sjPlot::tab_model(angle_model_u)
```
This shows that the angle model is the best one; only the direction of approach influences the chance of a U-turn. We also see that each of the three directions has a successively lower chance of U-turn. The odds ratio for ants coming from the main branch is one-third, for those coming from the primary branch it's even smaller, and for those coming from the secondary branch it's even smaller.

To interpret the model, let's look at the estimates for the coefficients for the fixed effect (which direction the ant is coming from). We can also use the model to make predictions for each possible approach.

```{r}
# use the table to make predictions from the model
approach.types$U.Turn.Prob <- predict(
  angle_model_u, newdata=approach.types, 
  re.form=NA, # this says ignore random effects in the predictions
  type="response" # this says convert output to probabilities
)
# approach.types <- approach.types %>%
#   mutate(Angle.From = case_when(
#            Angle.From == "main" ~ "Main",
#            Angle.From == "sharp" ~ "Secondary",
#            Angle.From == "shallow" ~ "Primary"))
approach.types %>% dplyr::select(Angle.From, Junction.Handedness, sharpIsLeft,
                          Sharp.Turn.Prob, U.Turn.Prob)
```

We can use this model to generate predictions for all turns, based just on the direction of approach. To do this, we need to join the table containing turns for all junctions (`turns`) with the table containing the predictions based on these factors (`approach.types`).

```{r}
turn_predictions <- left_join(turns, approach.types)
turn_predictions <- turn_predictions %>%
  mutate(Turn.Probability = case_when(isUTurn=="TRUE" ~
                                        U.Turn.Prob,
                                      isSharpTurn=="TRUE" ~
                                        (1-U.Turn.Prob)*Sharp.Turn.Prob,
                                      isSharpTurn=="FALSE" ~
                                        (1-U.Turn.Prob)*(1-Sharp.Turn.Prob)))
write.csv(turn_predictions, "data/NetworksSummer2021-TurnPredictions.csv",
          row.names=F)
```

Let's double check to make sure that the probabilities make sense: for every approach (the node it came from, and the node it's at), the sum of the turn probabilities (three different nodes it could go to) should be 1.

```{r}
turn_predictions %>%
  group_by(Node.From,Junction) %>%
  summarize(Total.Prob=sum(Turn.Probability))
```