R/buteo_age.R

In mottensmann/MyR: Collection of Custom Functions

#' Estimate the age of juvenile buzzards based on wing length
#'
#' @description Buzzard age is estimated based on a standard growth curve modeled from published data by Rob G. Bijlsma (1999)
#'
#' @param df data frame containing data
#' @param wing column name for wing measurements
#' @param sex column name giving the sex of individuals or NULL
#' @param unit unit of values. By default in mm
#'
#' @references
#' Bijlsma, RG 1999: Sex determination of nestling Common Buzzards Buteo buteo. Limosa 72, 1.10
#'
#' @examples
#' ## 240 mm wing length
#' buteo_age(data.frame(wing = c(240,240), sex = c("male","female")))
#'
#' @import ggplot2
#'
#' @export
#'
buteo_age <- function(df = NULL, wing = "wing", sex = "sex", unit = c("mm", "cm"), .plot = F) {
  unit <- match.arg(unit)

  ## get standard growth data
  path <- system.file("extdata", "buzzard_wing.txt", package = "MyR")
  data <- read.table(path, header = T, skip = 1,
                     colClasses = c("integer", "factor", "factor",
                                    "numeric", "numeric", "integer",
                                    "factor", "factor", "factor"))

  if (unit == "cm") {
    data[["mean"]] <- data[["mean"]]/10
    data[["sd"]] <- data[["sd"]]/10
  }

  ## build model
  if (is.null(sex)) {
    model <- lm(age ~ mean, data)
    fit.val <- data.frame(age = predict(model, data), mean = data[["mean"]])
  } else {
    model <- lm(age ~ mean + sex, data)
    fit.val <- data.frame(age = predict(model, data), mean = data[["mean"]], sex = data[["sex"]])
  }


  ## plot
  if (.plot == T) {
  plot <- ggplot(data, aes(y = age, x = mean, col = sex)) +
    geom_point(alpha = .8) +
    geom_line(data = fit.val, size = 1,
              aes(y = age, x = mean, col = sex)) +
    theme_classic(base_size = 12) +
    ylab("Age [days since hatching]") +
    xlab(paste0("Wing length [", unit, "]")) +
    scale_color_discrete(name = "",
                         breaks = c("1", "0"),
                         labels = c("Male", "Female")) +
    ggtitle("Common Buzzard growth curve",
            subtitle = "data: Bijlsma, RG (1999): Limosa 72") +
    annotate("text", y = 50, x = ifelse(unit == "mm", 70, 7), colour = "red",
             label = paste0("y=", round(coefficients(model)[[2]],2), "*x+",
                            round(coefficients(model)[[1]] - coefficients(model)[[3]],2),
                            " , R²=0.99")) +
    annotate("text", y = 47, x = ifelse(unit == "mm", 70, 7), colour = "blue",
             label = paste0("y=", round(coefficients(model)[[2]],2), "*x+",
                            round(coefficients(model)[[1]],2),
                            " , R²=0.99"))

  out <- plot
  } else {
    ## build model
    if (is.null(sex)) {
      model <- lm(age ~ mean, data)
      fit.val <- data.frame(age = predict(model, data), mean = data[["mean"]])
      ## format data
      to_predict <- data.frame(mean = df[[wing]])
    } else {
      model <- lm(age ~ mean + sex, data)
      fit.val <- data.frame(age = predict(model, data), mean = data[["mean"]], sex = data[["sex"]])
      ## format data
      to_predict <- data.frame(mean = df[[wing]],
                               sex = df[[sex]])
    }
    out <- predict.lm(model, to_predict, se.fit = T)
    return(out)
  }
  return(out)
}