R/PWD_get_gh.r

In ppwdeming: Precision Profile Weighted Deming Regression

#' Estimate of Variance Profile Functions (proportional)
#' @name PWD_get_gh
#'
#' @description
#' This code estimates the variance profiles, assumed proportional, of
#' the Rocke-Lorenzato form;
#' also provides the resulting weighted Deming fit and residuals.
#'
#' @usage
#' PWD_get_gh(X, Y, lambda=1, epsilon=1.e-8, printem=FALSE)
#'
#' @param X		the vector of predicate readings,
#' @param Y		the vector of test readings,
#' @param lambda		*optional* (default of 1) - the ratio of the X to the Y precision profile (defaults to 1),
#' @param epsilon		*optional* (default of 1.e-8) - convergence tolerance limit,
#' @param printem	  *optional* - if TRUE, routine will print out results as a `message`.
#'
#' @details
#' This workhorse routine optimizes the likelihood in the **unknown** *g*, *h*
#' setting over its *n*+4 parameters
#' (the two Rocke-Lorenzato precision profile parameters \eqn{\sigma}
#' and \eqn{\kappa}, the intercept \eqn{\alpha} and slope \eqn{\beta},
#' and the *n* latent true concentrations \eqn{\mu_i}).
#'
#' That is, the assumed forms are:
#'    * predicate precision profile model: \eqn{g_i = var(X_i) = \lambda\left(\sigma^2 + \left[\kappa\cdot \mu_i\right]^2\right)} and
#'    * test precision profile model: \eqn{h_i = var(Y_i) = \sigma^2 + \left[\kappa\cdot (\alpha + \beta\mu_i)\right]^2}.
#'
#' @returns A list containing the following components:
#'
#'   \item{alpha }{the fitted intercept}
#'   \item{beta }{the fitted slope}
#'   \item{fity }{the vector of predicted Y}
#'   \item{mu }{the vector of estimated latent true values}
#'   \item{resi }{the vector of residuals}
#'   \item{sigma }{the estimate of the Rocke-Lorenzato \eqn{\sigma}}
#'   \item{kappa }{the estimate of the Rocke-Lorenzato \eqn{\kappa}}
#'   \item{like }{the -2 log likelihood L}
#'
#' @author Douglas M. Hawkins, Jessica J. Kraker <krakerjj@uwec.edu>
#'
#' @examples
#' # library
#' library(ppwdeming)
#'
#' # parameter specifications
#' sigma <- 1
#' kappa <- 0.08
#' alpha <- 1
#' beta  <- 1.1
#' true  <- 8*10^((0:99)/99)
#' truey <- alpha+beta*true
#' # simulate single sample - set seed for reproducibility
#' set.seed(1039)
#' # specifications for predicate method
#' X     <- sigma*rnorm(100)+true *(1+kappa*rnorm(100))
#' # specifications for test method
#' Y     <- sigma*rnorm(100)+truey*(1+kappa*rnorm(100))
#'
#' # fit with RL precision profile to estimate parameters
#' RL_gh_fit  <- PWD_get_gh(X,Y,printem=TRUE)
#' # RL precision profile estimated parameters
#' cat("\nsigmahat=", signif(RL_gh_fit$sigma,6),
#'     "and kappahat=", signif(RL_gh_fit$kappa,6))
#'
#' @references Hawkins DM and Kraker JJ. Precision Profile Weighted Deming
#' Regression for Methods Comparison, on *Arxiv* (2025) <doi:10.48550/arXiv.2508.02888>
#'
#' @references Rocke DM, Lorenzato S (2012). A Two-Component Model for Measurement
#' Error in Analytical Chemistry.  *Technometrics*, **37:2**:176-184.
#'
#' @importFrom stats optimize
#'
#' @export

PWD_get_gh <- function(X, Y, lambda=1, epsilon=1.e-8, printem=FALSE) {

  key  <- order(X)
  sX   <- X[key]
  sY   <- Y[key]

  # Establish search ranges
  n    <- length(X)
  lowr <- 1:round(n/3)
  hir  <- round(2*n/3):n
  maxs <- max(abs(sY[lowr]-sX[lowr]))              # sigma
  maxk <- max(abs(log(sY[hir]/sX[hir])), na.rm=T)  # kappa
  gr   <- 1.618
  a    <- 0
  b    <- maxs
  fity <- X

  innr <- function(kappa, sigma, lambda) {
    do <- PWD_RL(X, Y, sigma, kappa, lambda)
    return(like=do$like)
  }

  h    <- maxs
  while (h > epsilon) {
    h        <- b - a
    c        <- b - h / gr
    d        <- a + h / gr
    vc       <- optimize(innr, c(0,maxk), c, lambda)
    fc       <- vc$objective
    vd       <- optimize(innr, c(0,maxk), d, lambda)
    fd       <- vd$objective
    if (fc < fd) {
      b <- d
    } else {
      a <- c
    }
  }
  sigma <- c
  kappa <- vc$minimum
  like    <- fc
  if (fd < fc) {
    sigma <- d
    kappa <- vd$minimum
    like    <- fd
  }
  do <- PWD_RL(X, Y, sigma, kappa, lambda)
  if (printem) {
    with(do, message(sprintf("alpha %6.3f beta %5.3f like %8.5g \n",
                             alpha, beta, like)))
  }
  return(list(alpha=do$alpha, beta=do$beta, fity=do$fity, mu=do$mu, resi=do$resi,
              sigma=sigma, kappa=kappa, like=do$like))
}