R/cal.wols.R

In bayclumpr: Bayesian Analysis of Clumped Isotope Datasets

#' Fit weighted OLS regression models on a given calibration dataset
#'
#' @param data The calibration dataset
#' @param replicates Number of bootstrap replicates
#' @param samples Number of samples per bootstrap replicate
#'
#' @return a \code{data.frame} with replicate-level
#'         regression parameters
#'
#' @importFrom stats lm
#'
#' @export

cal.wols <- function(data, replicates, samples = NULL) {
  if(is.null(samples)){samples = nrow(data)}
  reps <- lapply(1:replicates, function(x) {
    dataSub <- data[sample(seq_along(data[, 1]), nrow(data), replace = TRUE), ]
    Reg0 <- lm(D47 ~ Temperature, dataSub)
    wt <- 1 / lm(abs(Reg0$residuals) ~ Reg0$fitted.values)$fitted.values^2
    Reg <- summary(lm(D47 ~ Temperature, dataSub, weights = wt))
    res <- cbind.data.frame("alpha" = Reg$coefficients[1, 1], "beta" = Reg$coefficients[2, 1])
    return(res)
  })
  reps <- do.call(rbind, reps)
  return(reps)
}