# R/TDI.R In MethComp: Analysis of Agreement in Method Comparison Studies

#### Documented in TDI

#' Compute Lin's Total deviation index
#'
#' This index calculates a value such that a certain fraction of difference between methods will be numerically smaller than this. The TDI is a measure which esentially is a number K such that the interval [-K,K] contains the limits of agreement.
#'
#' If \code{boot==FALSE} a single number, the TDI is returned.
#' If \code{boot} is a number, the median and the 1-\code{alpha}/2 central interval
#' based on \code{boot} resamples are returned too, in a named vector of length 4.
#'
#' @param y1 Measurements by one method.
#' @param y2 Measurements by the other method.
#' @param p The fraction of items with differences numerically exceeding the TDI
#' @param boot If numerical, this is the number of bootstraps. If \code{FALSE} no confidence interval for the TDI is produced.
#' @param alpha 1 - confidende degree.
#'
#' @return A list with 3 components. The names of the list are preceeded by the
#'  criterion percentage, i.e. the percentage of the population that the TDI is
#'  devised to catch.
#'  \item{TDI}{The numerically computed value for the TDI. If \code{boot} is
#'             numeric, a vector of median and a bootstrap c.i. is appended.}
#'  \item{TDI}{The approximate value of the TDI}
#'  \item{Limits of Agreement}{Limits of agreement}
#' @references LI Lin: Total deviation index for measuring individual agreement with applications in laboratory performance and bioequivalence, Statistics in Medicine, 19, 255-270 (2000)
#' @author Bendix Carstensen, \email{bendix.carstensen@@regionh.dk}
#'
#' @examples
#'
#' data(plvol)
#' pw <- to.wide(plvol)
#'
#' @importFrom stats qnorm pnorm uniroot
#' @export
TDI <-
function( y1, y2, p = 0.05, boot=1000, alpha = 0.05 )
{
if( length(y1) != length(y2) )
stop( "Lengths of y1 and y2 must be the same!" )

# Analytical approximation
mu  <- mean(y1-y2)
sigma <- sd(y1-y2)
TDI.appr <- qnorm(1-p/2)*sqrt(mu^2+sigma^2)
names( TDI.appr ) <- "(approximate)"

# Limits of agreement
LoA <- mu + c(-1,1)*qnorm(1-p/2)*sigma
names( LoA ) <- c("lower","upper")

# Numerical calculation
FF <- function( x ) pnorm( (x-mu)/sigma ) - pnorm( (-x-mu)/sigma ) - (1-p)
int <- 2 * max(abs(LoA)) * c(-1,1)
TDI.num <- uniroot( FF, int )\$root
names( TDI.num ) <- "(numeric)"

# Bootstrap c.i.
if( is.numeric(boot) )
{
nn <- length(y1)
tdi <- numeric(boot)
for( i in 1:boot )
{
wh <- (1:nn)[sample(1:nn,nn,replace=T)]
tdi[i] <- TDI( y1[wh], y2[wh], p=p, boot=FALSE )[]
}
TDI.num <- c(TDI.num, quantile( tdi, c(0.5,alpha/2,1-alpha/2) ))
}

# Put it all together
res <- list( c(TDI.appr,TDI.num), LoA )
names( res ) <- c( paste((1-p)*100,"% TDI",sep=""),
paste((1-p)*100,"% Limits of Agreement",sep="") )
return( res )
}


## Try the MethComp package in your browser

Any scripts or data that you put into this service are public.

MethComp documentation built on Jan. 20, 2020, 1:12 a.m.