R/PlotVar.R

In gfaity/ReachStroke: Stroke reaching's Statistics

#' Plot Bland & Altman and Regression
#'
#' Plot the comparison of Ref and ToComp, with Ref as the reference
#' 1° left panel = B&A
#'     plot of Ref-ToComp as a function of (Ref+ToComp)/2
#' 2° rigth panel = regression of ToComp onto Ref (ToComp = a + bRef)
#'     plot of ToComp as a function of Ref
#'
#' @param Ref       Datas of reference system
#' @param ToComp    Datas of system to compare
#' @param U         "%"
#' @param RefT      Name of reference system
#' @param ToCompT   Name of system to compare
#' @param Cond      Main title
#' @param V         Current Variable
#' @param VT        Current Variable name
#' @param UnitCurr  Current Unit to use
#' @param Sym.Main
#'
#' @examples
#'  V <- Variables [iVariables]
#'  VT <- VariablesNames [iVariables]
#'  U <- "%"
#'  RefT <- RefNameToRead
#'  ToCompT <- ToCompNameToRead
#'  UnitCurr <- VariablesUnit[iVariables]
#'  Ref <- subset(D, App == RefName)[,V]
#'  ToComp <- subset(D, App == ToCompName)[,V]
#'  Cond <- paste(ToCompT, "VS", RefT)
#'
#' @export

PlotVar <- function(Ref,
                    ToComp,
                    U,
                    RefT,
                    ToCompT,
                    Cond,
                    V,
                    VT,
                    UnitCurr,
                    Sym.Main = "") {

  #################

  par.pty = par()$pty     # save current pty (to restore it in the end)
  par(pty = "s")          # make axis square

  #################

  if (UnitCurr  == "°"){
    Ref <- Ref*180/pi
    ToComp <- ToComp*180/pi
  }


  # Get the regression of ToComp onto Ref, ICC and B&A
  ToCompRef = lm(ToComp ~ Ref)                # plot(ToCompRef) for diagnostics
  sm = summary(ToCompRef)              # print(sm) for details

  slope = ToCompRef$coefficients[2]
  intercept = ToCompRef$coefficients[1]
  rsquared = sm$r.squared

  #library(irr)
  i = icc(cbind(Ref, ToComp))

  #library(BlandAltmanLeh)
  ba = bland.altman.stats(ToComp, Ref)

  #################

  # Display some information in the console
  cat(sprintf("\n--- %s : %s \n", Cond, V))
  cat(sprintf("Regression   :  %s = %0.2f*%s%+0.2f (r² = %0.2f)\n", ToCompT, slope, RefT, intercept, rsquared))
  cat(sprintf("%s : %s ICC = %.2f\n", i$type, i$model, i$value))
  cat(sprintf(
    "Bland & Altman : %s = %s%+1.2f,  SD = %1.2f (%s)\n",
    RefT,
    ToCompT,
    ba$mean.diffs,
    ba$critical.diff / ba$two, ## this is SD
    U
  ))

  #################

  # PLOT : Bland & Altman
  MinX = min((Ref + ToComp) / 2)
  MaxX = max((Ref + ToComp) / 2)
  DeltaY = (ba$upper.limit - ba$lower.limit) / 2
  Ylabel = sprintf("%s %s - %s %s", V, ToCompT, V, RefT)
  Xlabel = sprintf("%s mean", V, U)

  #if Variable name contains "Delta" then replace it by greek letter
  #Delta means difference End - Beg
  if (grepl("Delta", VT, fixed = TRUE)) {
    newVT <- gsub("Delta","", VT, fixed=TRUE)
    Ylabel = bquote(Delta ~ .(newVT) [.(ToCompT)] ~ .("-") ~ Delta ~ .(newVT) [.(RefT)] ~ .("(") * .(UnitCurr) * .(")"))
    Xlabel = bquote(Delta ~ .(newVT) ["mean"] ~ .("(") * .(UnitCurr) * .(")"))
    mainT = bquote(bold(Delta ~ .(paste(newVT, ":", Cond)) ))
  } else {
    Ylabel = bquote(.(VT) [.(ToCompT)] ~ .("-") ~ .(VT) [.(RefT)] ~ .("(") * .(UnitCurr) * .(")"))
    Xlabel = bquote(.(VT) ["mean"] ~ .("(") * .(UnitCurr) * .(")"))
    mainT = bquote(bold(.(paste(VT,":", Cond)) ))
  }

  bland.altman.plot(
    ToComp,
    Ref,
    xlim = c(MinX, MaxX),
    col = "blue",
    main = mainT,
    ylab = Ylabel,
    xlab = Xlabel
  )

  TheLabel = c(
    sprintf("%5.2f", ba$lower.limit),
    sprintf("%5.2f ±%0.2f", ba$mean.diffs, sd(ba$diffs)),
    sprintf("%5.2f", ba$upper.limit)
  )
  axis(
    side = 4,
    at = c(ba$lines),
    labels = TheLabel,
    las = 2
  )
  abline(0, 0, col = "gray", lty = 3)

  #################

  # PLOT : Regression of M onto Ref
  Xlabel = sprintf("%s from %s", V, RefT)
  Ylabel = sprintf("%s from %s", V, ToCompT)

  par(pty = "s")  # make axis square
  MaxXY = max(c(ToComp, Ref))
  MinXY = min(c(ToComp, Ref))

  #if Variable name contains "Delta" then replace it by greek letter
  if (grepl("Delta", VT, fixed = TRUE)) {
    newVT <- gsub("Delta","", VT, fixed=TRUE)
    Ylabel = bquote(Delta ~ .(newVT) [.(ToCompT)] ~ .("(") * .(UnitCurr) * .(")"))
    Xlabel = bquote(Delta ~ .(newVT) [.(RefT)] ~ .("(") * .(UnitCurr) * .(")"))
    mainT = bquote(bold(Delta ~ .(paste(newVT, ":", Cond)) ))
  } else {
    Ylabel = bquote(.(VT) [.(ToCompT)] ~ .("(") * .(UnitCurr) * .(")"))
    Xlabel = bquote(.(VT) [.(RefT)] ~ .("(") * .(UnitCurr) * .(")"))
    mainT = bquote(bold(.(paste(VT,":", Cond)) ))
  }

  plot(
    Ref,
    ToComp,
    col = "blue",
    main = mainT,
    xlab = Xlabel,
    ylab = Ylabel,
    xlim = c(MinXY, MaxXY),
    ylim = c(MinXY, MaxXY)
  )

  #Add ICC indication on plot
  text = paste("ICC =", i$type, i$model, round(i$value,2))
  mtext(text, side = 3, line = 0, outer = F, cex=0.7)


  # add the equation of the regression line (from )
  Equation = sprintf(
    "%s = %+0.2f %s %+0.2f\nr²=%0.2f",
    ToCompT,
    slope,
    RefT,
    intercept,
    rsquared
  )
  text(
    x = MinXY + (MaxXY - MinXY) * .04,
    y = MinXY + (MaxXY - MinXY) * .9,
    labels = Equation,
    adj = 0         # 0 = left ; 1 = right ; 0.5 = center
    #col = "blue"
  )
  abline(ToCompRef, col = "blue", lty = 2)  #regression line
  abline(c(0, 1), col = "black", lty = 3)   #identity line

  #################

  par(pty = par.pty)  # restore original pty

}