R/general_functions.R

Defines functions removeOutliers dctofm rc_age_err rcage d14c totErr intrErr sigma se rsd

Documented in d14c dctofm intrErr rcage rc_age_err removeOutliers rsd se sigma totErr

# General data reduction functions
# TODO: efficiency functions
# TODO: Calculate confidence interval
# TODO: better function for statistical normalization/standard (0 mean 1 var)
# TODO: Error propagation for D14C and rcage

#' Relative SD
#'
#' Provides the standard deviation of a vector normalized to the mean of the vector.
#'
#' @param x A numeric vector.
#'
#' @return The relative SD of x
#' @export
#'
rsd <- function(x) {
  sd(x)/mean(x)
}


#' Standard Error
#'
#' @param x A numeric vector.
#'
#' @return The SE of x
#' @export
#'
se <- function(x) {
  sqrt(var(x)/length(x))
}


#' Sigma
#'
#' Also known as a z-score. Provides an estimate of the quality of
#' a measurement relative to it's expected value.
#'
#' @param fmm Measured value.
#' @param fmc Expected value.
#' @param em Measurement error.
#' @param ec Error of expected value.
#'
#' @return A vector of sigma values
#' @export
#'
sigma <- function(fmm, fmc, em, ec = 0) {
  (fmm - fmc) / sqrt(em^2 + ec^2)
}

#' Normalized Value
#'
#' Normalize a value using the expected value or population mean.
#'
#' @param fmm Measured value.
#' @param fmc Expected value.
#'
#' @return A vector of normalized fm
#' @export
#'
normFm <- function (fmm, fmc) {
  #Calculate normalized fm
  (fmm - fmc) / fmc
}


#' Intrinsic Error
#'
#' Calculate the "extra" (intrinsic, residual) error.
#'
#' @param x A vector of measured values.
#' @param err A vector of measurement errors.
#' @param ... Additional parameters passed to sd.
#'
#' @return A vector of residual errors
#' @export
#'
intrErr <- function(x, err, ...) {
  stopifnot(is.numeric(x))
  stopifnot(is.numeric(err))
  xsd <- sd(x, ...)
  ifelse((xsd >= err),
         suppressWarnings(sqrt(xsd^2 - err^2)),
         NA)
}


#' Total error.
#'
#' Calculate the predicted variability in a population given a measurement
#' error and an intrisic error.
#'
#' @param rep_err Measurement error.
#' @param res_err Intrinsic error.
#' @param fm Fraction modern (optional).
#'
#' @return A vector of total errors.
#' @export
#'
totErr <- function(rep_err, res_err, fm) {
  stopifnot(is.numeric(rep_err))
  stopifnot(is.numeric(res_err))
  if (missing(fm)) {
    fm <- 1
  }
  stopifnot(is.numeric(fm))

  ifelse(is.na(res_err),
         rep_err,
         sqrt(rep_err^2 + (res_err * fm)^2)
  )
}

#' Calculate D14C.
#'
#' Calculate D14C given fraction modern and year of collection
#'
#'
#' @param fm Fraction modern.
#' @param yc Year of collection.
#'
#' @return A vector of D14C values
#' @export
#'
d14c <- function(fm, yc) {
  stopifnot(is.numeric(fm))
  stopifnot(is.numeric(yc))

  l <- 0.00012097 # 1/radiocarbon half life

  (fm * exp(l * (1950 - yc)) - 1) * 1000

}

#' Calculate radiocarbon years
#'
#' Calculate radiocarbon age from Fm
#'
#'
#' @param fm Fraction modern.
#'
#' @return radiocarbon age in years.
#' @export
#'
rcage <- function(fm) {
  stopifnot(is.numeric(fm))
  -8033 * log(fm)
}

#' Calculate radiocarbon age error
#'
#' @param fm Fraction modern
#' @param fm_err Error in the Fm
#'
#' @return radiocarbon age error in years
#' @export
#'
rc_age_err <- function(fm, fm_err) {
  sqrt((-8033 * 1/fm)^2 * (fm_err)^2)
}

#' Calculate Fm
#'
#' Calculate Fm from D14C
#'
#'
#' @param dc D14C of sample.
#' @param yc Year of collection.
#'
#' @return Fraction modern of sample.
#' @export
#'
dctofm <- function(dc, yc) {
  stopifnot(is.numeric(dc))
  stopifnot(is.numeric(yc))


  l <- 0.00012097 # 1/radiocarbon half life
  ((dc / 1000) + 1) / exp(l * (1950 - yc))

}

#' Remove outliers by interquartile range
#'
#' Outliers are defined as points more than a multiplier * IQR above the
#' third quartile or multiplier * IQR below the first quartile. The
#' default multiplier is 1.5.
#'
#' @param x A vector of numeric values
#' @param multiplier A value to multiply by IQR
#'
#' @return The vector with outliers replaced by `NA`
#' @export
#'
removeOutliers = function(x, multiplier = 1.5) {
  stopifnot(is.numeric(x),
            is.numeric(multiplier),
            length(multiplier) == 1)
  # Get Q1 and Q3
  qnt <-  quantile(x, probs=c(.25, .75))

  iqt <- multiplier * IQR(x)

  # Apply on a copy of the original data
  y = x
  y[x < (qnt[1] - iqt)] = NA
  y[x > (qnt[2] + iqt)] = NA
  y
  # Remove incomplete cases and return the resulted variable
  #y[complete.cases(y)]
}
blongworth/amsdata documentation built on June 15, 2022, 5:35 p.m.