library(ggplot2)


my_theme<-function(base_size = 14) {
  theme_bw(base_size = base_size,
           base_family = "serif") %+replace%
    theme(
      # The whole figure
      plot.title = element_text(size = rel(1), face = "bold",
      margin = margin(0,0,5,0), hjust = 0.5),

      # Area where the graph is located
      panel.grid.minor = element_blank(),
      panel.border = element_blank(),

      # The axes
      axis.title = element_text(size = rel(0.85), face = "bold"),
      axis.text = element_text(size = rel(0.70), face = "bold"),
      axis.line = element_line(color = "black"),
  
      # The legend
      legend.title = element_text(size = rel(0.85), face = "bold"),
      legend.text = element_text(size = rel(0.70), face = "bold"),
      legend.key = element_rect(fill = "transparent", colour = NA),
      legend.key.size = unit(1.5, "lines"),
      legend.background = element_rect(fill = "transparent", colour = NA),

      # The labels in the case of facetting
      strip.background = element_rect(fill = "#17252D", color = "#17252D"),
      strip.text = element_text(size = rel(0.85), face = "bold",
      color = "white", margin = margin(5,0,5,0))
    )
}

# Change the default theme
theme_set(my_theme())
# package names
packages<-c("tidyverse", "here", "mcp", "lubridate", "knitr", "ezknitr", "loo", "flextable")

# install any packages not previously installed
installed_packages<-packages %in% rownames(installed.packages())
if(any(installed_packages == FALSE)){
  install.packages(packages[!installed_packages])
}

# load packages
invisible(lapply(packages, library, character.only = TRUE))

Load observation data

obs<-read_csv(here("DRUM_materials/!revised_materials/ODBA_nesting/data/raw_data/all_0A_obs.csv"))
head(obs)
## # A tibble: 6 x 14
##    ...1 ID    datetime  time   status  behavior  cygnets ncygnets adults nadults
##   <dbl> <chr> <chr>     <time> <chr>   <chr>     <chr>      <dbl> <chr>  <chr>  
## 1     1 0A    4/1/2021  13:20  On nest Loafing ~ No            NA Marke~ <NA>   
## 2     2 0A    4/25/2021 01:40  On nest Incubati~ No            NA Marke~ <NA>   
## 3     3 0A    4/27/2021 10:54  On nest Incubati~ No            NA Marke~ <NA>   
## 4     4 0A    4/29/2021 09:43  On nest Incubati~ No            NA Marke~ <NA>   
## 5     5 0A    5/2/2021  13:45  On nest Incubati~ No            NA Marke~ <NA>   
## 6     6 0A    5/4/2021  11:19  On nest Incubati~ No            NA Marke~ <NA>   
## # ... with 4 more variables: comments <chr>, foraging <dbl>, loafing <dbl>,
## #   incubating <dbl>

Load raw accelerometry data

df<-read_csv(here("DRUM_materials/!revised_materials/ODBA_nesting/data/raw_data/0A_sensor_data.csv"))
head(df)
## # A tibble: 6 x 23
##   device_id UTC_datetime        UTC_date   UTC_time datatype satcount U_bat_mV
##       <dbl> <dttm>              <date>     <time>   <chr>       <dbl>    <dbl>
## 1    201692 2021-04-01 00:00:10 2021-04-01 00:00:10 GPS            10     4183
## 2    201692 2021-04-01 00:14:41 2021-04-01 00:14:41 GPS            11     4186
## 3    201692 2021-04-01 00:29:41 2021-04-01 00:29:41 GPS            11     4183
## 4    201692 2021-04-01 00:44:41 2021-04-01 00:44:41 GPS            10     4180
## 5    201692 2021-04-01 00:59:41 2021-04-01 00:59:41 GPS             6     4180
## 6    201692 2021-04-01 01:14:49 2021-04-01 01:14:49 GPS             6     4180
## # ... with 16 more variables: bat_soc_pct <dbl>, solar_I_mA <dbl>, hdop <dbl>,
## #   Latitude <dbl>, Longitude <dbl>, Altitude_m <dbl>, speed_km_h <dbl>,
## #   direction_deg <dbl>, temperature_C <dbl>, mag_x <dbl>, mag_y <dbl>,
## #   mag_z <dbl>, acc_x <dbl>, acc_y <dbl>, acc_z <dbl>, ...23 <lgl>

Deal with bursts.

The current transmitter settings are GPS fixes every 15 minutes throughout the 24 hour cycle and ACC fixes every 5 minutes with a 10Hz 3 second burst.

Convert to POSIXct format

df$UTC_datetime<-ymd_hms(df$UTC_datetime)

Filter out data before ACC sensors were activated and from the start of observations of nesting

df<-df %>% filter(UTC_datetime>"2021-04-24")

Assign a burst ID to all rows

df<-df %>% 
  group_by(five_min_burst=cut(UTC_datetime, "5 min")) %>% 
  mutate(burst_row_number=1:n(),
         burst_ID=cur_group_id())

Convert acc units into g-force

df$acc_x_g<-df$acc_x/1024
df$acc_y_g<-df$acc_y/1024
df$acc_z_g<-df$acc_z/1024

Calculate ODBA

df<-df %>% 
  group_by(burst_ID) %>% 
  mutate(ODBA=sum(abs(acc_x_g-mean(acc_x_g)),abs(acc_y_g-mean(acc_y_g)),abs(acc_z_g-mean(acc_z_g))))

df$julian_day<-yday(df$UTC_datetime)

Average hourly ODBA

df<-df %>% 
  mutate(hourly_burst=cut(UTC_datetime, "1 hour")) %>% 
  group_by(hourly_burst) %>% 
  mutate(hr_odba=mean(ODBA))

df$time<-ymd_hms(df$hourly_burst)

Subset to just 1 value per hour

hourly<-df %>% 
  distinct(hourly_burst, hr_odba)

Join nesting status observations

obs$Day<-as.Date(obs$datetime, format="%m/%d/%Y")
obs<-obs %>% 
  mutate(date_time=ymd_hms(paste(Day,time)))
obs$Jday<-yday(obs$date_time)
obs$hour<-hour(obs$date_time)

hourly$date<-ymd_hms(hourly$hourly_burst)
hourly$Jday<-yday(hourly$date)
hourly$hour<-hour(hourly$date)


hourly<-left_join(x=hourly, y=obs, by=c("Jday", "hour"))

The marked swan was consistently seen incubating the next from April 25 - May 13. On May 25th and May 28th the collared swan was not seen, but was presumed to be on the nest. June 2nd was the first date that cygnets were observed. I’ll take the midpoint between those two times as the estimate of hatch date from human observers.

Plot both datasets with aligned x-axes.

hourly$hours_since_start<-as.numeric(difftime(hourly$date, hourly$date[1], units="hours"))

raw1<-ggplot(hourly, aes(hours_since_start, hr_odba)) + geom_point() + 
  labs(x="Hourly time intervals", y="Average hourly ODBA value", title="Accelerometry Dataset")+
  my_theme()

obs$new_date<-mdy_hms(paste(obs$datetime, obs$time, sep=" "))
first_cyg<-head(obs[obs$status=="Cygnets","new_date"],1)
last_nest<-tail(obs[obs$status=="On nest","new_date"],1)

hatch<-as.POSIXct((as.numeric(first_cyg$new_date)+
                     as.numeric(last_nest$new_date))/2,
                  origin="1970-01-01")

raw2<-ggplot(hourly, aes(date, ID)) +
  geom_point(data=filter(hourly, is.na(status) !=TRUE), aes(date, col = status), size =4)+
  geom_vline(aes(xintercept = hatch),linetype="dashed", size=1.5, color="purple")+
  labs(x="Date", color="Nesting status", title="Visual Observations Dataset")+
   scale_color_manual(labels=c("Cygnets present","Incubating"), values=c("blue", "orange"))+
  my_theme()+
  theme(legend.position = c(0.1,0.9))

cowplot::plot_grid(raw1, raw2, align = "h", axis="l", ncol=1)

Use mcp to detect the breakpoint for the time of nest hatching

Models can take some time to fit, so the following 2 chunks won’t be run in this .Rmd, and instead I’ll import a saved model object.

hourly<-drop_na(hourly, "hours_since_start")

m0<-list(
  hr_odba~1)

m1<-list(
 hr_odba~1,
        ~1
)

m2<-list(
  hr_odba~1,
         ~1,
         ~1
)

f0<-mcp(m0,hourly[,c("hr_odba", "hours_since_start")],
         par_x="hours_since_start")

f1<-mcp(m1, hourly[,c("hr_odba", "hours_since_start")],
         par_x="hours_since_start")

f2<-mcp(m2,hourly[,c("hr_odba", "hours_since_start")],
         par_x="hours_since_start")
mods<-list(f0, f1, f2)
for(i in 1:length(mods)){
  mods[[i]]$loo<-loo(mods[[i]])
}

loo_compare(mods[[1]]$loo,
            mods[[2]]$loo,
            mods[[3]]$loo)
mods<-readRDS(here("DRUM_Materials/!revised_materials/ODBA_nesting/saved_models/all_with_loo.rds"))

loos<-loo_compare(mods[[1]]$loo, mods[[2]]$loo, mods[[3]]$loo)
loo_df<-as.data.frame(loos[1:3, 1:2])
loo_df<-cbind.data.frame('Model Syntax'=c("Two Intercepts", 
                                          "Three Intercepts",
                                          "One Intercept"), loo_df)


loo_df %>%  flextable() %>% 
  colformat_double(digits=1) %>% 
  set_header_labels(values=list("Model Syntax" = "Model Syntax", 
                         "elpd_diff" = "ELPD Difference", 
                        "se_diff" = "SE of Difference")) %>% 
  set_caption(caption = "Leave-one-out Cross Validation was used to compute the Estimated Log Predictive Density (ELPD) of three different models; one fit with a single intercept, one with two intercepts, and one with three intercepts. The higher ELPD for the model with two intercepts indicates that it has a higher predictive accuracy.") %>% 
  border_outer() %>% 
  vline(part="header") %>% 
  align(align ='center', part = 'all') %>% 
  bg(part="header", bg="#dbdbdb") %>% 
  bold(part="header") %>%  
  set_table_properties(layout="autofit")

Top model has a single breakpoint and two segments with separate intercepts for activity levels.

summary(mods[[2]])
## Family: gaussian(link = 'identity')
## Iterations: 9000 from 3 chains.
## Segments:
##   1: hr_odba ~ 1
##   2: hr_odba ~ 1 ~ 1
## 
## Population-level parameters:
##     name  mean lower upper Rhat n.eff
##     cp_1 893.4 891.0 895.8    1  3928
##    int_1   3.9   3.7   4.2    1  5956
##    int_2   7.5   7.3   7.7    1  5382
##  sigma_1   3.6   3.5   3.7    1  5691
plot_pars(mods[[2]])

In the the figure below, the vertical purple dashed line is the estimated hatch time based on visual observations, and the vertical brown solid line is the estimated hatch time based on the piecewise regression model fit. Therefore the piecewise regression model is able to detect the breakpoint at the time of nest hatch when activity transitions from a lower state to a higher state as incubation ends.

invisible({capture.output({out<-summary(mods[[2]])})})
mu<-out$mean[[1]]

invisible({capture.output({fitted<-as.data.frame(fitted(mods[[2]]))})})
names(fitted)[[2]]<-"hr_odba"
fitted$date<-as.POSIXct(fitted$hours_since_start*60*60, 
                        origin = "2021-04-24 05:00:00")
ggplot(hourly, aes(date, hr_odba)) + geom_point() + 
  labs(x="Date", y="Average hourly ODBA value", color="Nesting Status")+
  geom_line(data=fitted,aes(date, hr_odba), color="grey", size=2)+
  geom_errorbar(data=fitted, aes(ymin=Q2.5, ymax=Q97.5), color="red", 
                linetype="dashed", size=0.4)+
  geom_rug(data=filter(hourly, is.na(status) !=TRUE), 
           aes(date, hr_odba, col = status), 
           size =2, sides="b")+
  scale_color_manual(labels=c("Cygnets present","Incubating"), values=c("blue", "orange"))+
  geom_vline(aes(xintercept = hatch), linetype="dashed", size=1.5, color="purple")+
    geom_vline(aes(xintercept = as.POSIXct(mu*3600, 
                                  origin = "2021-04-24 05:00:00")),
                   linetype="solid", color="brown", size=1.5)+
  my_theme()+
  theme(legend.position = c(0.1, 0.9))+

  annotate("segment", x=as.POSIXct(mu*3600, origin = "2021-04-24 05:00:00")+10000,
           xend= as.POSIXct(mu*3600, origin = "2021-04-24 05:00:00")+750000, y=20,
           yend=22.5 , color="brown", size=2)+
  
  annotate(x=as.POSIXct(mu*3600, origin = "2021-04-24 05:00:00")+750000,
           y=+Inf, label="Piecewise Regression", vjust=2, geom="label", color="brown")+
   
  annotate("segment", x=hatch-500000, xend=hatch-10000, y=22.5, yend=20, color="purple", size=2)+
  annotate(x=hatch-750000, y=+Inf, label="Visual Observations", vjust=2, geom="label", color="purple")
knitr::include_graphics(here("writing/methods_paper/images/nesting2.jpg"))

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] flextable_0.6.8 loo_2.4.1       ezknitr_0.6     knitr_1.30     
##  [5] lubridate_1.8.0 mcp_0.3.2       here_0.1        forcats_0.5.1  
##  [9] stringr_1.4.0   dplyr_1.0.8     purrr_0.3.4     readr_2.1.2    
## [13] tidyr_1.2.0     tibble_3.1.6    tidyverse_1.3.1 ggplot2_3.3.5  
## 
## loaded via a namespace (and not attached):
##  [1] matrixStats_0.57.0   fs_1.5.0             bit64_4.0.5         
##  [4] httr_1.4.2           rprojroot_1.3-2      tensorA_0.36.2      
##  [7] tools_4.0.2          backports_1.1.10     bslib_0.3.1         
## [10] utf8_1.1.4           R6_2.4.1             ggdist_3.1.1        
## [13] DBI_1.1.1            colorspace_1.4-1     withr_2.3.0         
## [16] tidyselect_1.1.2     bit_4.0.4            compiler_4.0.2      
## [19] cli_3.2.0            rvest_1.0.2          HDInterval_0.2.2    
## [22] arrayhelpers_1.1-0   xml2_1.3.2           officer_0.4.0       
## [25] labeling_0.3         posterior_1.2.1      sass_0.4.0          
## [28] scales_1.1.1         checkmate_2.0.0      ggridges_0.5.3      
## [31] systemfonts_0.2.3    digest_0.6.28        rmarkdown_2.11      
## [34] base64enc_0.1-3      pkgconfig_2.0.3      htmltools_0.5.2     
## [37] dbplyr_2.1.1         fastmap_1.1.0        rlang_1.0.2         
## [40] readxl_1.3.1         rstudioapi_0.13      farver_2.0.3        
## [43] svUnit_1.0.3         jquerylib_0.1.4      generics_0.0.2      
## [46] jsonlite_1.8.0       vroom_1.5.7          distributional_0.3.0
## [49] zip_2.1.1            magrittr_2.0.1       bayesplot_1.8.0     
## [52] patchwork_1.1.1      Rcpp_1.0.8           munsell_0.5.0       
## [55] fansi_0.4.1          abind_1.4-5          gdtools_0.2.2       
## [58] lifecycle_1.0.1      stringi_1.5.3        yaml_2.2.1          
## [61] plyr_1.8.6           grid_4.0.2           parallel_4.0.2      
## [64] crayon_1.5.0         lattice_0.20-41      haven_2.3.1         
## [67] hms_1.1.1            pillar_1.7.0         uuid_0.1-4          
## [70] reshape2_1.4.4       codetools_0.2-16     reprex_2.0.1        
## [73] glue_1.6.2           evaluate_0.14        data.table_1.13.0   
## [76] modelr_0.1.8         vctrs_0.3.8          tzdb_0.2.0          
## [79] cellranger_1.1.0     gtable_0.3.0         assertthat_0.2.1    
## [82] xfun_0.26            broom_0.7.12         coda_0.19-3         
## [85] tidybayes_3.0.2      ellipsis_0.3.2