# R/BA.plot.r In MethComp: Analysis of Agreement in Method Comparison Studies

#### Documented in BA.plot

```#' Bland-Altman plot of differences versus averages.
#'
#' For two vectors of equal length representing measurements of the same
#' quantity by two different methods, the differences are plotted versus the
#' average. The limits of agreement (prediction limits for the differences) are
#' plotted, optionally a regression of differences of means is given too. Works
#'
#' A plot of the relationship between the methods is produced; either a
#' Bland-Altman plot of the differences versus averages, or a 45 degree
#' rotation as a conversion between the methods. If \code{model=NULL} a simple
#' regression of averages on differences is made by calling \code{DA.reg}, and
#' the specified conversion plotted.
#'
#' @param y1 Numerical vector of measurements by 1st method. Can also be a
#' @param y2 Numerical vector of measurements by 2nd method. Must of same
#' length as \code{x}. Ignored if a \code{\link{Meth}} or a
#' \code{\link{MethComp}} objects is given for \code{y1}.
#' @param meth.names Label for the method names.
#' @param wh.comp Which methods should be compared. Either numerical or
#' character.
#' @param pl.type What type of plot should be made, \code{"BA"} for differences
#' versus averages, \code{"conv"} for method 1 versus method 2.
#' @param dif.type How should difference depend on the averages. \code{"const"}
#' or \code{"lin"}.
#' @param sd.type How should the standard deviation depend on the averages.
#' \code{"const"} or \code{"lin"}.
#' @param model Should a variance component model be used to compute the limits
#' of agreement? If \code{NULL} a simple analysis is made; other possibilities
#' @param eqax Should the axes be identical? If a Bland-Altman plot is drawn,
#' the axis for the differences will have the same extent as the axis for the
#' averages, but centered on 0 (see \code{diflim}).
#' @param axlim The limits of the axes.
#' @param diflim The limits of the difference axis.
#' @param grid Should a grid be drawn? If numeric it indicates the places where
#' the grid is drawn.
#' @param N.grid How many grid-lines should be drawn.
#' @param col.grid What should be the color of the grid?
#' @param points Logical. Should the observed points be drawn?
#' @param col.points What color should they have?
#' @param cex.points How large should they be?
#' @param pch.points What plot character for the points
#' @param lwd Numerical vector of 3, giving the width of the conversion line
#' (mean difference) and the limits of agreement.
#' @param col.lines What color should the lines have.
#' @param repl.conn Should replicate measurements be connected (within items)?
#' @param col.conn Color of connecting lines.
#' @param lwd.conn Width of connecting lines.
#' @param xlab x-axis label.
#' @param ylab y-axis label.
#' @param eqn Logical. Should the equations linking the methods be shown on the
#' plot? If a Bland-Altman plot is made, both the equations linking the methods
#' and the equation for the differences versus the averages are shown.
#' @param col.eqn Color for equations
#' @param font.eqn Font for equations
#' @param digits How many digits after the decimal point should be used when
#' showing the equations.
#' @param Transform Transformation applied to data prior to analysis. Plots are
#' made on the original scale after back-transformation.
#' @param mult Logical. If TRUE, ratios of measurement instead of differences
#' will be plotted in the Bland-Altman plot on a logarithmic axis, and limits
#' of agreement will be given on this scale?  This gives the same analysis as
#' using \code{Transform="log"}, but a different plot. Using another
#' transformation than the log is accommodated, but no LoA is shown on the
#' axis.
#' @param alpha 1 minus the confidence level. If \code{NULL} a multiplier of 2
#' is used for constructing prediction limits, otherwise a t-quantile with d.f.
#' equal th number of items minus 1.
#' @param ... Further parameters passed on to \code{\link{plot.MethComp}}
#' @return An object of class \code{\link{MethComp}} and either \code{DA.reg}
#' (if \code{model=NULL}) or \code{BA.est} (if \code{model} is character).
#' @author Bendix Carstensen \email{bendix.carstensen@@regionh.dk},
#' \url{http://BendixCarstensen.com}.
#' @references JM Bland and DG Altman: Statistical methods for assessing
#' agreement between two methods of clinical measurement, Lancet, i, 1986, pp.
#' 307-310.
#'
#' JM Bland and DG Altman. Measuring agreement in method comparison studies.
#' Statistical Methods in Medical Research, 8:136-160, 1999.
#'
#' B Carstensen: Comparing methods of measurement: Extending the LoA by
#' regression. Stat Med. 2010 Feb 10;29(3):401-10.
#' @keywords models design
#' @examples
#'
#' data( ox )
#' ox <- Meth( ox )
#' # The simplest possible Bland-Altman plot
#' BA.plot( ox )
#'
#' ## With bells and whistles, comparing the naive and model
#' par( mfrow=c(2,2) )
#' BA.plot( ox, model=NULL, repl.conn=TRUE, col.lines="blue",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=TRUE, dif.type="lin", pl.type="BA", sd.type="lin",
#'          grid=1:9*10, digits=3,font.eqn=1)
#' par(new=TRUE)
#' BA.plot( ox, model="linked", repl.conn=TRUE, col.lines="red",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=FALSE, dif.type="lin", pl.type="BA", sd.type="lin",
#'         grid=1:0*10, digits=3)
#' BA.plot( ox, model=NULL, repl.conn=TRUE, col.lines="blue",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=TRUE, dif.type="lin", pl.type="conv", sd.type="lin",
#'         grid=1:9*10, digits=3,font.eqn=1)
#' par(new=TRUE)
#' BA.plot( ox, model="linked", repl.conn=TRUE, col.lines="red",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=FALSE, dif.type="lin", pl.type="conv", sd.type="lin",
#'          grid=1:9*10, digits=3)
#' # The same again, but now logit-transformed
#' BA.plot( ox, model=NULL, repl.conn=TRUE, col.lines="blue",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=TRUE, dif.type="lin", pl.type="BA", sd.type="lin",
#'          grid=1:9*10, digits=3,font.eqn=1,Transform="pctlogit")
#' par(new=TRUE)
#' BA.plot( ox, model="linked", repl.conn=TRUE, col.lines="red",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=FALSE, dif.type="lin", pl.type="BA", sd.type="lin",
#'          grid=1:0*10, digits=3,Transform="pctlogit")
#' BA.plot( ox, model=NULL, repl.conn=TRUE, col.lines="blue",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=TRUE, dif.type="lin", pl.type="conv", sd.type="lin",
#'          grid=1:9*10, digits=3,font.eqn=1,Transform="pctlogit")
#' par(new=TRUE)
#' BA.plot( ox, model="linked", repl.conn=TRUE, col.lines="red",
#'          axlim=c(0,100), diflim=c(-50,50), xaxs="i", yaxs="i",
#'          las=1, eqn=FALSE, dif.type="lin", pl.type="conv", sd.type="lin",
#'          grid=1:9*10, digits=3,Transform="pctlogit")
#'
#' @rdname plot
#' @export
BA.plot <-
function( y1, y2, meth.names = NULL,
wh.comp = 1:2,
pl.type = "BA",
dif.type = "const",
sd.type = "const",
model = if( inherits(y1,"Meth") & has.repl(y1) ) "exch"
else NULL,
eqax = FALSE,
axlim = if( is.data.frame(y1) ) range(y1\$y) else range(c(y1,y2)),
diflim = NULL,
grid = TRUE,
N.grid = 10,
col.grid = grey(0.9),
points = TRUE,
col.points = "black",
cex.points = 1,
pch.points = 16,
lwd = c(3,1,1),
col.lines = "blue",
repl.conn = FALSE,
col.conn = "gray",
lwd.conn = 1,
xlab = NULL,
ylab = NULL,
eqn = FALSE,
col.eqn = col.lines,
font.eqn = 2,
digits = 2,
Transform = if( mult ) "log" else NULL,
mult = FALSE,
alpha = NULL,
... )
{
# Allow sloppy definition of arguments
pl.type <- ifelse( tolower( substr( pl.type,1,1) ) == "c", "conv" , "BA"  )
dif.type <- ifelse( tolower( substr(dif.type,1,1) ) == "c", "const", "lin" )
sd.type <- ifelse( tolower( substr( sd.type,1,1) ) == "c", "const", "lin" )
if( !is.null(model) )
model <- ifelse( tolower( substr(   model,1,1) ) == "l", "linked", "exch" )

if( is.vector( y1 ) )
# If we have a vector, check if it has the same length as the second argument
if( length(y1)!=length(y2) )
stop( "Arguments y1 and y2 must have same length, but",
"length(y1)=",  length(y1),
", length(y2)=",  length(y2) )
else
{
# And if they are of same length, make a Meth object out of it using
# supplied names if givem
tmp <- data.frame(y1,y2)
if( is.character(meth.names) ) names(tmp) <- meth.names
y1 <- Meth( tmp, y=1:2, print=FALSE )
}

if( is.data.frame( y1 ) )
{
# If the dataframe is not a Meth object make it, Meth will take
# care of the possible errors if the right columns are not there.
if( !inherits(y1,"Meth") )
y1 <- Meth( y1, print=FALSE )
# Select the two methods to compare and subset the Meth object to
# the two methods that we plot, and make sure wh.comp hold their names
if( is.numeric(wh.comp) ) wh.comp <- levels(y1\$meth)[wh.comp]
obj <- Meth( y1[y1\$meth %in% wh.comp,], print=FALSE )
}

else stop("Wrong data structrue for y1 supplied: str(y1):", str(y1) )

# So we turn this into a MethComp object
if( is.null(model) )
M.obj <- DA.reg( obj, DA.slope = dif.type=="lin",
Transform = Transform )
else
Transform = Transform )

# Then we compute various parameters for the plotting

# axlim:
if( is.null(axlim) ) axlim <- range( obj\$y )

# diflim:
if( is.null(diflim) ) diflim <- axlim - mean(axlim)

# And then we can use the default machinery to plot this:
plot.MethComp( M.obj,
wh.comp = wh.comp,
pl.type = pl.type,
dif.type = dif.type,
sd.type = sd.type,
axlim = axlim,
diflim = diflim,
points = points,
grid = grid,
N.grid = N.grid,
col.grid = col.grid,
lwd = lwd,
col.lines = col.lines,
col.points = col.points,
pch.points = pch.points,
repl.conn = repl.conn,
col.conn = col.conn,
lwd.conn = lwd.conn,
eqn = eqn,
col.eqn = col.eqn,
font.eqn = font.eqn,
digits = digits,
alpha = alpha,
mult = mult,
... )

attr( M.obj, "pl.type" ) <- c( pl.type =  pl.type,
dif.type = dif.type,
sd.type =  sd.type,
model = model )
invisible( M.obj )
}

######################################################################
# Utility to connect points to means.

connect2mean <-
function( obj, wh.comp,
pl.type = "conv",
col.conn = "gray",
lwd.conn = 1,
... )
{
wob <- to.wide( obj )
# The points
wob <- cbind( wob[,wh.comp[2]],
wob[,wh.comp[1]],
# - and the item averages
ave(wob[,wh.comp[2]],wob[,"item"]),
ave(wob[,wh.comp[1]],wob[,"item"]) )
# Convert to D-A coordinates if required
if( pl.type == "BA" )
{
c(0.5, 1))