### LoA.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: apr 5 2019 (09:46)
## Version:
## Last-Updated: sep 7 2021 (11:35)
## By: Brice Ozenne
## Update #: 102
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * LOA - documentation
##' @title Limits of Agreement
##' @description Computes the limits of agreement between two measurements.
##' @name LoA
##'
##' @param x [numeric vector] vector containing the first measurement, each value corresponding to an individual.
##' @param y [numeric vector] vector containing the second measurement, each value corresponding to an individual (same as in \code{x}).
##' @param conf.level.CI [numeric 0-1] confidence level of the confidence intervals.
##' @param conf.level.LoA [numeric 0-1] confidence level for the limits of agreements.
##'
##' @details The bias is estimated as the average difference between y and x, i.e. bias = E[y-x].
##' The upper/lower limit of the limit of agreement are estimated as the bias plus/minus the LoA.multiplier times the standard deviation of y-x,
##' i.e. E[y-x] +/- LoA.multiplier * sqrt(Var[y-x]).
##'
##' The standard errors, confidence intervals, and p-values are computed for the limit of the limit of agreement are computed using a delta method.
##'
##'
##'
##'
## * LoA - example
##' @examples
##'
##' set.seed(10)
##' X <- rnorm(100)
##' Y <- X + rnorm(100)
##'
##' e.loa <- LoA(x = X, y = Y)
##' e.loa
##' autoplot(e.loa)
##'
##' if(require(blandr)){ ## same results as blandr.statistics
##' e.blandr <- blandr.statistics(method1 = Y, method2 = X)
##' e.loa$LoA$estimate - c(e.blandr$bias, e.blandr$lowerLOA, e.blandr$upperLOA)
##' e.loa$LoA$lower - c(e.blandr$biasLowerCI, e.blandr$lowerLOA_lowerCI, e.blandr$upperLOA_lowerCI)
##' e.loa$LoA$upper - c(e.blandr$biasUpperCI, e.blandr$lowerLOA_upperCI, e.blandr$upperLOA_upperCI)
##' }
## * LoA - code
##' @rdname LoA
##' @export
LoA <- function(x, y, conf.level.CI = 0.95, conf.level.LoA = 0.95){
## ** check arguments
if(length(x) != length(y)){
stop("Argument \'x\' and \'y\' must have same lenght \n")
}
## ** normalize arguments
alpha <- 1 - conf.level.CI
LoA.multiplier <- qnorm(conf.level.LoA)
n <- length(x)
df <- n-1
critical.quantile <- qt(1-alpha/2, df = df)
data <- data.frame(x = x,
y = y,
difference = y - x,
average = (y+x)/2)
data.LoA <- data.frame(matrix(NA, nrow = 4, ncol = 5,
dimnames = list(c("bias","lowerLoA","upperLoA","halfWidthLoA"),
c("estimate","se","lower","upper","p.value"))))
## ** compute quantities of interest
mean.diff <- mean(data$difference)
var.diff <- stats::var(data$difference)
## var(sigma^2) = 2*sigma^4/(n-1)
var.var.diff <- 2*var.diff^2/(n-1)
sd.diff <- sqrt(var.diff)
## var(sigma) = var(sqrt(sigma^2)) = var(sigma^2)/(4*sigma^2) by a delta method
var.sd.diff <- var.var.diff / (4*var.diff)
se.diff <- sd.diff/sqrt(n)
## var(mu + q * sigma) = var(mu) + q^2 * var(sigma) + 2 * 0 (no covariance between mu and sigma)
var.LoA <- (se.diff)^2 + LoA.multiplier^2 * var.sd.diff
se.LoA <- sqrt(var.LoA)
## bias
data.LoA["bias","estimate"] <- mean.diff
data.LoA["bias","se"] <- se.diff
## lower LoA
data.LoA["lowerLoA","estimate"] <- mean.diff - LoA.multiplier * sd.diff
data.LoA["lowerLoA","se"] <- se.LoA
## upper LoA
data.LoA["upperLoA","estimate"] <- mean.diff + LoA.multiplier * sd.diff
data.LoA["upperLoA","se"] <- se.LoA
## half width
data.LoA["halfWidthLoA","estimate"] <- LoA.multiplier * sd.diff
data.LoA["halfWidthLoA","se"] <- LoA.multiplier * sqrt(var.sd.diff)
## p-value and CI
data.LoA[,"lower"] <- data.LoA[,"estimate"] - critical.quantile * data.LoA[,"se"]
data.LoA[,"upper"] <- data.LoA[,"estimate"] + critical.quantile * data.LoA[,"se"]
data.LoA[,"p.value"] <- 2*(1-pt(abs(data.LoA[,"estimate"]/data.LoA[,"se"]), df = df))
data.LoA["halfWidthLoA","p.value"] <- NA
## ** export
out <- list(data = data,
LoA = data.LoA,
conf.level.CI = conf.level.CI,
conf.level.LoA = conf.level.LoA)
class(out) <- "LoA"
return(out)
}
## * print.LoA
print.LoA <- function(x, digits = 2, eps = 1e-4, print = TRUE, ...){
out <- x$LoA
out$estimate <- round(out$estimate, digits = digits)
out$se <- round(out$se, digits = digits)
out$lower <- round(out$lower, digits = digits)
out$upper <- round(out$upper, digits = digits)
out$p.value <- format.pval(out$p.value, digits = digits, eps = eps)
out$lower <- paste0("[",out$lower, ";",out$upper,"]")
out$upper <- NULL
names(out)[names(out)=="lower"] <- "CI"
if(print == TRUE){
print(out)
}
return(invisible(out))
}
## * autoplot.LoA
##' @title Bland and Altman Plot with Limits of Agreement
##' @description Display a Bland and Altman plot with limits of agreement.
##'
##' @param object An object of class \code{LoA}, output of the \code{LoA} function.
##' @param display.ci [logical] should the confidence intervals be displayed?
##' @param plot [logical] should the plot be printed?
##' @param line.size [numeric >0] the width of the horizontal lines.
##' @param alpha [numeric 0-1] transparency parameter for the confidence intervals.
##' @param name.legend [character] character string displayed in the caption.
##' @param colors [character of length 3] colors used to display the bias, lower, and upper LoA.
##' @param labels [character of length 3] text shown in the caption for the bias, lower, and upper LoA.
##' @param linetype [numeric of length 3] type of lines used to display the bias, lower, and upper LoA.
##' @param ... not used, for compatibility with the generic method.
##'
##' @export
autoplot.LoA <- function(object, display.ci = TRUE, plot = TRUE,
line.size = 1.5, alpha = 0.25, name.legend = "",
labels = c("bias","LoA (lower)","LoA (upper)"),
colors = c("red","blue","forestgreen"), linetype = c(1,2,2),
...){
if(display.ci){
keep.type <- c("estimate","lower","upper")
}else{
keep.type <- c("estimate")
}
object$LoA$type <- rownames(object$LoA)
## display
gg <- ggplot2::ggplot()
if(display.ci){
LoA <- object$LoA[c("bias","lowerLoA","upperLoA"),]
x.range <- range(object$data$average)
gg <- gg + ggplot2::geom_rect(data = LoA, aes_string(ymin = "lower", ymax = "upper", fill = "type"), alpha = alpha, xmin = x.range[1], xmax = x.range[2])
}else{
LoA <- object$LoA[c("bias"),]
}
gg <- gg + ggplot2::geom_point(data = object$data, mapping = ggplot2::aes_string(x = "average", y = "difference"))
gg <- gg + ggplot2::geom_hline(data = data.frame(type = LoA$type, estimate = LoA$estimate),
ggplot2::aes_string(yintercept = "estimate", linetype = "type", color = "type"), size = line.size)
gg <- gg + ggplot2::scale_colour_manual(name = name.legend,
labels = labels,
values = c("bias" = colors[1],
"lowerLoA" = colors[2],
"upperLoA" = colors[3]))
gg <- gg + ggplot2::scale_fill_manual(name = name.legend,
labels = labels,
values = c("bias" = colors[1],
"lowerLoA" = colors[2],
"upperLoA" = colors[3]))
gg <- gg + ggplot2::scale_linetype_manual(name = name.legend,
labels = labels,
values = c("bias" = linetype[1],
"lowerLoA" = linetype[2],
"upperLoA" = linetype[3]))
if(plot == TRUE){
print(gg)
}
## export
return(invisible(gg))
}
######################################################################
### LoA.R ends here
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.