A function to calculate and plot Jacobson-Truax based reliable clinical change. Here as a gist rather than a package because it has not been fully checked yet (especially for crit B and C). Input & comments welcome!

JT_RCC.R

JT_RCC <- function(x.pre, x.post, norm.M = NA, norm.SD = NA, RC.crit ="auto", reliability, ppid, plot = T, plottype = "JT", higherIsBetter = F) {
  
  SEmeasurement <- sd(x.pre) * sqrt(1-reliability) # standard error of measurement based on the sd of the pre measurement/baseline 
  
  Sdiff<- sqrt(2*SEmeasurement^2) # the SD of the SEm for a difference score
  
  # determine the criterion to determine recovery:
  
  ### C is preferred, followed by B, then A.. 
  
  if(is.na(norm.M)|is.na(norm.SD))  { # if no norm data is available, only crit A is  possible
    crittype <- "crit A"
    crit <- mean(x.pre) - 1.96 * sd(x.pre)  # #The level of functioning at post should fall outside the range of the clinical population (more than 1.96 standard deviation in the 'more healthy' direction). 
  } 
  else { # if norm data is availabe:
    
    if (norm.M + 3* norm.SD < mean(x.pre) - 3* sd(x.pre)) { # if the pre/clinical sample does not overlap with the norm sample, use criterion C
      crittype <- "crit C"
      crit <- ( (sd(x.pre) * norm.M) + (norm.SD * mean(x.pre)) ) /(sd(x.pre) + norm.SD) #The level of functioning at post should place the patient closer to the mean of the comparison norm data than the mean of the clinical norm data / pre measurement - weigted by SD of both clin and norm "
    } 
    
    else { # if the distributions of the pre/clinical sample overlaps with the norm distribution, use criterion B:
      crittype <- "crit B"
      crit <- norm.M + 1.96 * norm.SD # The level of functioning at post should fall within the range of the comparison non-clinical group,i.e. within 1.96 standard deviation of the mean of the comparison norm data  
    }
  }
  
  # see -> http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.615.8373&rep=rep1&type=pdf
  
  # time to compute things: 
  JTdf <- NULL
  JTdf$ppid <- ppid
  x.pre <- as.numeric(x.pre)
  x.post <- as.numeric(x.post)
  JTdf$pre <- x.pre
  JTdf$post <- x.post
  JTdf <- as.data.frame(JTdf)
  
  if(higherIsBetter) {JTdf$change <- JTdf$pre - JTdf$post } else
  {JTdf$change <- JTdf$post - JTdf$pre}
  
  # write values to JTdf:
  JTdf$SEmeasurement <- SEmeasurement
  JTdf$Sdiff<- Sdiff 
  JTdf$crittype <- crittype
  JTdf$crit <- crit
  
  # compute Reliable Change Index:
  JTdf$RCI <- JTdf$change / Sdiff
  
  # compute Jacobson-Truax classification:
  JTdf$JTclass_RCI <- NA
  JTdf$JTclass_RCI [JTdf$post <= crit & JTdf$RCI <= -1.96 ] <- "recovered"
  JTdf$JTclass_RCI [JTdf$post <= crit & JTdf$RCI > -1.96 ]  <- "non reliably recovered"
  JTdf$JTclass_RCI [JTdf$post > crit & JTdf$RCI <= -1.96 ]  <- "improved"
  JTdf$JTclass_RCI [JTdf$post > crit & JTdf$RCI > -1.96 ]   <- "unchanged"
  JTdf$JTclass_RCI [JTdf$RCI >= 1.96 ]                      <- "deteriorated"
  
  JTdf$JTclass_RCI <- factor(JTdf$JTclass_RCI, levels = rev(c("recovered", "non reliably recovered", "improved", "unchanged", "deteriorated")))
  
  cat("Jacobson-Truax classification: \n")
  print(table(JTdf$JTclass_RCI))
  
  # previous version also computed Wise table 1 (https://pdfs.semanticscholar.org/5e34/1ae1a311602280bebd31bfe98b5d17ed3ca1.pdf) based classifications but I've cut them out for the moment.  
  #
  JTdf <<- JTdf
  
  if(plot) {
    if (plottype == "JT") {
      JTdf$JTclassPlot <- JTdf$JTclass_RCI
      
      for (l in 1: length(levels(JTdf$JTclassPlot))) {
        levels(JTdf$JTclassPlot)[l] <- paste0(levels(JTdf$JTclassPlot)[l], ": ", table(JTdf$JTclassPlot)[[l]])
      }
      
      # created vector of plot colours linked to specific level names:  
      cols <- c("orchid3", "olivedrab3", "skyblue3", "coral1", "palevioletred3")
      
      names(cols) <- c(levels(JTdf$JTclassPlot)[grep("^recovered", levels(JTdf$JTclassPlot))], levels(JTdf$JTclassPlot)[grep("improved", levels(JTdf$JTclassPlot))], levels(JTdf$JTclassPlot)[grep("unchanged", levels(JTdf$JTclassPlot))], levels(JTdf$JTclassPlot)[grep("deteriorated",  levels(JTdf$JTclassPlot))], levels(JTdf$JTclassPlot)[grep("non reliably", levels(JTdf$JTclassPlot))]) 
      #
      
      if(sum(is.na(JTdf$post > 0))) {  
        warning( sum(is.na(JTdf$post)), " cases have missing post data - these are ommitted from the plot") }
      
      
      require(ggplot2)
      ggplot(JTdf[!(is.na(JTdf$post > 0)),], aes(x=pre, y=post, colour = JTclassPlot)) + 
        geom_point( size = 2.5) + 
        geom_abline(intercept = 0, slope = 1) + 
        geom_abline(intercept = sqrt(2*SEmeasurement^2) * 1.96 , slope = 1, linetype =3) + 
        geom_abline(intercept = sqrt(2*SEmeasurement^2) * -1.96 , slope = 1, linetype =3) +
        geom_hline(yintercept = crit, linetype = 2) +
        theme_classic() + 
        theme(panel.grid.major = element_line(color = "gray95")) + 
        labs(title = "Jacobson-Truax plot", x = "pre", y = "post", colour = "Jacobson-Truax \n Classification:") + 
        scale_colour_manual(values = cols)
      
    }
  }
  
}