# Electronic Supplementary Material to 
# Lendvai et al. Analysis of the optimal duration of behavioral observations based on an 
# automated continuous monitoring system in tree swallows (Tachycineta bicolor): is one 
# hour good enough?



# read in nest visit data
load("Visitdata_ESM.RData")

# extracting whole day feeding rates
ffv.day  <-  visitdata$Fday #female
names(ffv.day)  <- row.names(visitdata)
mfv.day  <-  visitdata$Mday #male 
names(mfv.day)  <- row.names(visitdata)

visitdata$Fday  <- visitdata$Mday  <-  NULL

make.long = function(data, firsthour = 7L, numhours = 15L, numtime = 4L) {
  # reshaping data into a long format
  tmp1  <- data
  names(tmp1) = substr(names(tmp1), 2, 6)
  numcol  <- numhours*numtime
  tmp2  <-   rbind(tmp1[,1:numcol], tmp1[,(numcol + 1):(2*numcol)])
  tmp2$nestID  <- factor(row.names(tmp1))
  tmp2$sex  <- factor(rep(c("F", "M"), each = nrow(tmp1)))
  
  tmp3  <-  reshape(tmp2, direction = "long", idvar = c("nestID", "sex"),  varying = 1:numcol, sep = "/")
  tmp4 <- reshape(tmp3, direction = "long", idvar = c("nestID", "sex", "time"), 
                  varying = list(paste0(firsthour:(firsthour+numhours-1), "h")), 
                  timevar = "starthour", times = firsthour:(firsthour+numhours-1), 
                  v.names = "visit")
  
  visitlong  <-  tmp4[order(tmp4$nestID, tmp4$starthour, tmp4$time),]
}

visitlong  <-  make.long(data = visitdata)


# function to calculate R2 based on subsamples and the whole-day feeding rate
# hours: at which hour should the subsample be started? given as an integer
# howlong: how long should the subsample be? Given as an integer in 15 min increments,
#           therefore for a 1h-sample, give 4
# the function assumes the structure and variable names of the dataframe provided here ('visitlong')


getr2  <-  function(hours, howlong, data, plot = T, ...) 
{
  r <- list(fN = integer(0), mN = integer(0),fr2 = numeric(0), mr2 = numeric(0))
  
  for (h in hours) {
    maxtime <- max(data$time)  
    extrahour <- howlong %/% maxtime
    extratime <- howlong %% maxtime
    tmp1 <- data[data$starthour >= h & data$starthour <= (h + extrahour),]
    tmp2 <- tmp1[tmp1$starthour < (h + extrahour) | tmp1$time <= extratime, ]
    fsample <- tmp2[tmp2$sex == "F",]
    msample  <- tmp2[tmp2$sex == "M",]
    
    ffv.sample  <-  tapply(fsample$visit, fsample$nestID, sum, na.rm = F)/
      tapply(fsample$visit, fsample$nestID, length)*maxtime
    mfv.sample  <-  tapply(msample$visit, msample$nestID, sum)/
      tapply(msample$visit, msample$nestID, length)*maxtime  
    
    ffv.day2  <-  ffv.day[match(names(ffv.sample), names(ffv.day))]
    mfv.day2  <-  mfv.day[match(names(mfv.sample), names(mfv.day))]
    
    # R-sqared
    fr2  <-  summary(lm(ffv.day2~ffv.sample))$r.squared
    mr2  <-  summary(lm(mfv.day2~mfv.sample))$r.squared
    
    i  <-  h-(min(hours)-1)
    r$fN[[i]]  <-  length(ffv.sample[!is.na(ffv.sample)]) 
    r$mN[[i]]  <-  length(ffv.sample[!is.na(mfv.sample)]) 
    r$fr2[[i]]  <-  fr2 ; r$mr2[[i]]  <-  mr2
  }
  
  if (plot == T) {  
    fy  <-  r$fr2; my  <-  r$mr2 ; ylab  <-  expression(paste("", R^2))
    plot(hours, fy, type = "l", col = "red", ylim = c(0, 1), 
         xlab = "Start of sampling", ylab = ylab) 
    lines(hours, my, col = "blue")
  }
return(r)
}

# 1h-samples throughout the day
getr2(hours = 7:21, howlong = 4, data = visitlong)


# repeat it with varying observation durations
subsamples <- lapply(7:17, function(hrs) {
  lapply(1:16, function(x) getr2(hours = hrs, howlong = x, data = visitlong, plot = F))
})

# create objects to keep the results in
# females
fr2.matrix <- sapply(subsamples, function(hrs) sapply(hrs, function(obslength) obslength$fr2))
# males
mr2.matrix <- sapply(subsamples, function(hrs) sapply(hrs, function(obslength) obslength$mr2))

r2.df <- data.frame(fr2 = as.vector(t(fr2.matrix)), mr2 = as.vector(t(mr2.matrix)),
                   obslength = rep(seq(0.25, 4, 0.25), each = length(subsamples)))

# Plot the R2 in function of the observation duration. Each line corresponds to a starting hour
# females
matplot(fr2.matrix, type = "l", ylab = expression(paste("", R^2)), 
        xlab = "Duration of behavioral sampling (h)", 
        xaxt = "n", col = "red", lty = 1, ylim = c(0,1))
axis(1, at = seq(0, dim(fr2.matrix)[1], 4), labels = 0:4)

# males
matplot(mr2.matrix, type = "l", ylab = expression(paste("", R^2)), 
        xlab = "Duration of behavioral sampling (h)", 
        xaxt = "n", col = "blue", lty = 1, ylim = c(0,1))
axis(1, at = seq(0, dim(fr2.matrix)[1], 4), labels = 0:4)

# optimization function suggested by Froud & Abel 2014
# doi:10.1371/journal.pone.0114468
giveqss.optimum  <- function(curve = NULL, obslength = NULL, r2 = NULL) {
  qss  <- function(x, y) (1-x)^2 + (1-y)^2 
  if (is.null(obslength)) {
    x <- unique(curve$origData$obslength)
    y <- unique(predict(curve))
  }
  else if (is.null(curve)) {
    x <- unique(obslength)
    y <- unique(round(r2, 6))
  }
  else stop ('specify either the curve object or both obslength and r2!')
  opt  <- qss(x, y) 
  # it gives the same result, if we rescale the original variable as in Froud & Abel:
  # qss((((curve$origData$obslength/4) -1)*-1), predict(curve)))
  list(x = x[which.min(opt)], y = y[which.min(opt)])
} 

library(drc)
fcurve.MM <- drm(fr2~obslength, data = r2.df, fct = MM.3(names = c("lower", "upper", "mid")))
plot(fcurve.MM, log = "")
summary(fcurve.MM)

mcurve.MM <- drm(mr2~obslength, data = r2.df, fct = MM.3(names = c("lower", "upper", "mid")))
plot(mcurve.MM, log = "")
summary(mcurve.MM)

# optimum value for females
(fopt.MM  <- giveqss.optimum(fcurve.MM))

# optimum value for males
(mopt.MM  <- giveqss.optimum(mcurve.MM))

# plot it
op <- par(mar = c(5,5,2,2))
plot(predict(fcurve.MM)~r2.df$obslength, xlim = c(0,4), ylim = c(0.0, 0.8), 
     type = "l", col = "red", 
     xlab = "Duration of behavioral sampling (h)", 
     ylab = expression(R^2))
lines(r2.df$obslength, predict(mcurve.MM), col = "blue")
legend("bottomright", col = c("red", "blue"), lty = 1, legend = c("female", "male"))
# Add minor tick marks
require(Hmisc)
minor.tick(nx=4)

# add to the plot above
points(fopt.MM, col = "red", pch = 19)
points(mopt.MM, col = "blue", pch = 19)
arrows(x0 = fopt.MM$x, y0 = fopt.MM$y, x1 = fopt.MM$x, y1 = -0.025, col = "red", lty = "dashed", length = 0.2)
arrows(x0 = mopt.MM$x, y0 = mopt.MM$y, x1 = mopt.MM$x, y1 = -0.025, col = "purple", lty = "dashed", length = 0.2)