#Loading Packages

library(ggplot2)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(report)
library(ggpubr)

#Reading in the Data

dovedata <- read.csv("Bird Data - Sheet1.csv")

head(dovedata)
##   N NumBird time_s Dist Pecking Peck_rate
## 1 1       2     NA  5.8      NA        NA
## 2 2       2     NA 30.7      NA        NA
## 3 3       2     NA  9.5      NA        NA
## 4 4       1     60 16.2       2 0.1234568
## 5 5       1     NA 44.8      NA        NA
## 6 6       2     NA 40.1      NA        NA
str(dovedata)
## 'data.frame':    27 obs. of  6 variables:
##  $ N        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ NumBird  : int  2 2 2 1 1 2 1 3 8 4 ...
##  $ time_s   : num  NA NA NA 60 NA NA NA NA NA NA ...
##  $ Dist     : num  5.8 30.7 9.5 16.2 44.8 40.1 30.2 19.3 34.4 20.5 ...
##  $ Pecking  : int  NA NA NA 2 NA NA NA NA NA NA ...
##  $ Peck_rate: num  NA NA NA 0.123 NA ...

#Setting GGplot theme

Asia_Theme <- theme_classic()+
              theme(axis.text.x = element_text(size = 16, hjust =.75),
                    plot.title = element_text(size=16, face='bold'),
  axis.text.y = element_text(size = 16),
  axis.title.x = element_text(size = 18),
  axis.title.y = element_text(size = 18))

#Model 1: Escape Distance by Group Size

#Looking at the effect of group size on escape distance
mod.dist <- lm(Dist ~ NumBird, data = dovedata)
summary(mod.dist)
## 
## Call:
## lm(formula = Dist ~ NumBird, data = dovedata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.4936  -4.6674   0.2232   5.6183  20.2183 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24.8698     2.9236   8.507 7.53e-09 ***
## NumBird      -0.2881     0.8155  -0.353    0.727    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.205 on 25 degrees of freedom
## Multiple R-squared:  0.004968,   Adjusted R-squared:  -0.03483 
## F-statistic: 0.1248 on 1 and 25 DF,  p-value: 0.7268
plot(mod.dist)

# Plotting the Data

plot1 <- ggplot(data = dovedata, aes(x= NumBird, y = Dist))+
  geom_point(size=3)+
  theme_classic()+
  labs(
       x= "Group size",
       y= "Flight Initiation Distance")+
  Asia_Theme

plot1

#Model 2: Pecking Rate by Group Size

#Model of pecking rate by group size
mod.pecks <- lm(Peck_rate ~ NumBird, data = dovedata)
summary(mod.pecks)
## 
## Call:
## lm(formula = Peck_rate ~ NumBird, data = dovedata)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.092510 -0.027785 -0.003786  0.032736  0.139778 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.16246    0.06545   2.482   0.0421 *
## NumBird     -0.02307    0.02723  -0.847   0.4248  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07485 on 7 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.09302,    Adjusted R-squared:  -0.03655 
## F-statistic: 0.7179 on 1 and 7 DF,  p-value: 0.4248

Plotting the Data

plot2 <- ggplot(data = dovedata, aes(x= NumBird, y = Peck_rate))+
  geom_point(size = 3)+
  labs(
       x= "Group size",
       y= "Peck Rate")+
    theme_classic()+
  xlim(0.5, 4.5)+
  Asia_Theme

plot2
## Warning: Removed 18 rows containing missing values (`geom_point()`).

Simulated Data

set.seed(1)
Beta.1 <- -1
fake.grp <- runif(20, min = 1, max = 10)
error <-rnorm(20, mean = 1, sd = 2)
y1 <- Beta.1*fake.grp + error

simulate.data <- data.frame(x = fake.grp, y = y1)
simulate.plot <- ggplot(simulate.data, aes(x = fake.grp, y = y1))+
  geom_point(size = 3)+
  geom_smooth(method = lm, se = FALSE)+
  labs(title = "Prediction",
       x = "Group Size",
       y = "Flight Initiation Distance")+
   Asia_Theme+
    theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank()  #remove y axis ticks
        )

  
simulate.plot
## `geom_smooth()` using formula = 'y ~ x'

set.seed(1)
Beta.1 <- -1
fake.grp <- runif(20, min = 1, max = 10)
error <-rnorm(20, mean = 1, sd = 2)
y1 <- fake.grp + error

simulate.data <- data.frame(x = fake.grp, y = y1)
simulate.plot <- ggplot(simulate.data, aes(x = fake.grp, y = y1))+
  geom_point(size = 3)+
  geom_smooth(method = lm, se = FALSE)+
    labs(title = "Prediction",
       x = "Group Size",
       y = "Peck Rate")+
  Asia_Theme+
      theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank()  #remove y axis ticks
        )

  
simulate.plot
## `geom_smooth()` using formula = 'y ~ x'

Regrouping data

dovedata$group[dovedata$NumBird == 1] <- "individual"
dovedata$group[dovedata$NumBird == 2] <- "pair"
dovedata$group[dovedata$NumBird >= 3] <- "group"

dovedata
##     N NumBird time_s Dist Pecking  Peck_rate      group
## 1   1       2     NA  5.8      NA         NA       pair
## 2   2       2     NA 30.7      NA         NA       pair
## 3   3       2     NA  9.5      NA         NA       pair
## 4   4       1  60.00 16.2       2 0.12345679 individual
## 5   5       1     NA 44.8      NA         NA individual
## 6   6       2     NA 40.1      NA         NA       pair
## 7   7       1     NA 30.2      NA         NA individual
## 8   8       3     NA 19.3      NA         NA      group
## 9   9       8     NA 34.4      NA         NA      group
## 10 10       4     NA 20.5      NA         NA      group
## 11 11       7     NA 20.2      NA         NA      group
## 12 12       9     NA 22.5      NA         NA      group
## 13 13       1     NA 10.7      NA         NA individual
## 14 14       5     NA 24.2      NA         NA      group
## 15 15       2  60.00 23.5       2 0.03333333       pair
## 16 16       4  60.25 20.4       4 0.06639004      group
## 17 17       2  42.00 25.5       1 0.02380952       pair
## 18 18       1  29.00 34.4       0         NA individual
## 19 19       1  24.00 31.9       0         NA individual
## 20 20       2  82.00 17.8      21 0.25609756       pair
## 21 21       1  30.20 30.2       0         NA individual
## 22 22       5  43.00 18.8       0         NA      group
## 23 23       1  27.00 27.3       5 0.18518519 individual
## 24 24       2  49.38 28.1       0         NA       pair
## 25 25       3  30.63 20.4       3 0.09794319      group
## 26 26       2  22.59 15.2       2 0.08853475       pair
## 27 27       3  31.75 26.7       4 0.12598425      group

Looking at the effect of group size on escape distance with categorical predictors

mod.dist.2 <- lm(Dist ~ group, data = dovedata)
summary(mod.dist.2)
## 
## Call:
## lm(formula = Dist ~ group, data = dovedata)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17.51  -3.69  -0.24   3.83  18.30 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.740      2.834   8.024 2.99e-08 ***
## groupindividual    5.473      4.251   1.287    0.210    
## grouppair         -0.940      4.118  -0.228    0.821    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.962 on 24 degrees of freedom
## Multiple R-squared:  0.09457,    Adjusted R-squared:  0.01912 
## F-statistic: 1.253 on 2 and 24 DF,  p-value: 0.3036
mod.dist.3 <- lm(Peck_rate ~ group, data = dovedata)
summary(mod.dist.3)
## 
## Call:
## lm(formula = Peck_rate ~ group, data = dovedata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.07663 -0.03086 -0.01191  0.02921  0.15565 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     0.096772   0.046206   2.094   0.0811 .
## groupindividual 0.057548   0.073058   0.788   0.4608  
## grouppair       0.003671   0.061125   0.060   0.9541  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08003 on 6 degrees of freedom
##   (18 observations deleted due to missingness)
## Multiple R-squared:  0.1112, Adjusted R-squared:  -0.1851 
## F-statistic: 0.3752 on 2 and 6 DF,  p-value: 0.7022
plot3 <-dovedata %>%
  mutate(group = fct_relevel(group, 
            "individual", "pair", "group")) %>%
  ggplot(aes(x= group, y = Dist))+
  geom_boxplot()+
  labs(
       x= "Group type",
       y= "Flight Initiation Distance")+
  Asia_Theme

plot3

plot4 <-  dovedata %>%
  mutate(group = fct_relevel(group, 
            "individual", "pair", "group")) %>%
  ggplot(aes(x= group, y = Peck_rate))+
  geom_boxplot()+
  labs(x= "Group type",
       y= "Peck Rate")+
    Asia_Theme

plot4
## Warning: Removed 18 rows containing non-finite values (`stat_boxplot()`).

#Compiling Plots into one

panelfig <-ggarrange(plot1,plot2,plot3,plot4)
## Warning: Removed 18 rows containing missing values (`geom_point()`).
## Warning: Removed 18 rows containing non-finite values (`stat_boxplot()`).
panelfig

ggsave("plot1_doves.png", panelfig, width = 10, height = 7, units = "in")

Histogram of group sizes

hist.plot <-ggplot(data = dovedata, aes(x = NumBird))+
  geom_histogram(bins = 9)
hist.plot+
      labs(y = "Observations",
       x = "Group Size")+
  scale_x_continuous(breaks=c(1,2,3,4,5,6,7,8,9))+
  Asia_Theme