R/cdc.R

In scan: Single-Case Data Analyses for Single and Multiple Baseline Designs

#' Conservative Dual-Criterion Method
#'
#' The `cdc()` function applies the Conservative Dual-Criterion Method (Fisher,
#' Kelley, & Lomas, 2003) to scdf objects. It compares phase B data points to
#' both phase A mean and trend (OLS, bi-split, tri-split) with an additional
#' increase/decrease of .25 SD. A binomial test against a 50/50 distribution is
#' computed and p-values below .05 are labeled "systematic change".
#'
#' @inheritParams .inheritParams
#' @param trend_method Method used to calculate the trend line. Default is
#'   `trend_method = "OLS"`. Possible values are: `"OLS"`, `"bisplit"`, and
#'   `"trisplit"`. `"bisplit"`, and `"trisplit"` should only be used for cases
#'   with at least five data-points in both relevant phases.
#' @param conservative The CDC method adjusts the original mean and trend lines
#'   by adding (expected increase) or subtracting (expected decrease) an
#'   additional .25 SD before evaluating phase B data. Default is the CDC method
#'   with `conservative = .25`. To apply the Dual-Criterion (DC) method, set
#'   `conservative = 0`.
#' @return \item{cdc}{CDC Evaluation based on a p-value below .05.}
#'   \item{cdc_exc}{Number of phase B datapoints indicating expected change.}
#' \item{cdc_nb}{Number of phase B datapoints.} \item{cdc_p}{P value of Binomial
#' Test.} \item{cdc_all}{Overall CDC Evaluation based on all instances/cases of
#' a Multiple Baseline Design.} \item{N}{Number of cases.}
#'  \item{decreasing}{Logical argument from function call (see `Arguments`
#'   above).}
#' \item{conservative}{Numeric argument from function call (see `Arguments`
#' above).} \item{case_names}{Assigned name of single-case.} \item{phases}{-}
#' @author Timo Lueke
#' @references Fisher, W. W., Kelley, M. E., & Lomas, J. E. (2003). Visual Aids
#'   and Structured Criteria for Improving Visual Inspection and Interpretation
#'   of Single-Case Designs. *Journal of Applied Behavior Analysis, 36*,
#'   387-406. https://doi.org/10.1901/jaba.2003.36-387
#' @family overlap functions
#' @keywords overlap
#' @examples
#'
#' ## Apply the CDC method (standard OLS line)
#' design <- design(n = 1, slope = 0.2)
#' dat <- random_scdf(design, seed = 42)
#' cdc(dat)
#'
#' ## Apply the CDC with Koenig's bi-split and an expected decrease in phase B.
#' cdc(exampleAB_decreasing, decreasing = TRUE, trend_method = "bisplit")
#'
#' ## Apply the CDC with Tukey's tri-split, comparing the first and fourth phase.
#' cdc(exampleABAB, trend_method = "trisplit", phases = c(1,4))
#'
#' ## Apply the Dual-Criterion (DC) method (i.e., mean and trend without
#' ##shifting).
#' cdc(exampleAB_decreasing, decreasing = TRUE, trend_method = "bisplit", conservative = 0)
#'
#'
#' @export
cdc <- function(data, 
                dvar, 
                pvar, 
                mvar, 
                decreasing = FALSE, 
                trend_method = c("OLS", "bisplit", "trisplit"), 
                conservative = .25, 
                phases = c(1, 2)) {

  
  
  
  check_args(
    by_class(decreasing, "logical"),
    by_call(trend_method, "cdc"),
    within(conservative, 0, 1)
  )
  
  trend_method <- trend_method[1]
  
  # set attributes to arguments else set to defaults of scdf
  if (missing(dvar)) dvar <- dv(data) else dv(data) <- dvar
  if (missing(pvar)) pvar <- phase(data) else phase(data) <- pvar
  if (missing(mvar)) mvar <- mt(data) else mt(data) <- mvar
  
  data  <- .prepare_scdf(data, na.rm = TRUE)
  data  <- recombine_phases(data, phases = phases)$data
  
  N       <- length(data)
  cdc_na  <- rep(NA, N)  # total data points in phase A
  cdc_nb  <- rep(NA, N)  # total data points in phase B
  cdc_exc <- rep(NA, N)  # exceeding data points in phase B
  cdc     <- rep(NA, N)  # CDC rule evaluation of change
  cdc_p   <- rep(NA, N)  # binomial p (50/50)
  cdc_all <- NA          # CDC rule evaluation of all "cases"
  
  for(i in 1:N) {
    A <- data[[i]][data[[i]][, pvar] == "A",]
    B <- data[[i]][data[[i]][, pvar] == "B",]
    cdc_na[i] <- nrow(A)
    cdc_nb[i] <- nrow(B)

    if ((cdc_na[i] < 5 || cdc_nb[i] < 5) && trend_method != "OLS") {
      stop(
        "The selected method for trend estimation should not be applied ",
        "with less than five data points per phase.\n"
      )
    }

    if(trend_method == "bisplit"){
      x <- A[[mvar]]
      y <- A[[dvar]]
      # na.rm = FALSE for now to prevent misuse; will draw no line if NA present
      md1 <- c(
        (median(y[1:floor(length(y) / 2)], na.rm = FALSE)),
        median(x[1:floor(length(x) / 2)], na.rm = FALSE)
      )
      md2 <- c(
        (median(y[ceiling(length(y) / 2 + 1):length(y)], na.rm = FALSE)),
        median(x[ceiling(length(x) / 2 + 1):length(x)], na.rm = FALSE)
      )
      md <- as.data.frame(rbind(md1, md2))
      names(md) <- c(dvar,mvar)
      formula <- as.formula(paste0(dvar, "~", mvar))
      model <- lm(formula, data = md, na.action = na.omit)
    }
    
    if(trend_method == "trisplit"){
      x <- A[[mvar]]
      y <- A[[dvar]]
      # na.rm = FALSE for now to prevent misuse; will draw no line if NA present
      md1 <- c(
        (median(y[1:floor(length(y) / 3)], na.rm = FALSE)),
        median(x[1:floor(length(x) / 3)], na.rm = FALSE)
      )
      md2 <- c(
        (median(y[ceiling(length(y) / 3 * 2 + 1):length(y)], na.rm = FALSE)),
        median(x[ceiling(length(x) / 3 * 2 + 1):length(x)], na.rm = FALSE)
      )
      md <- as.data.frame(rbind(md1, md2))
      names(md) <- c(dvar,mvar)
      formula <- as.formula(paste0(dvar,"~",mvar))
      model <- lm(formula, data = md, na.action = na.omit)
    }
    
    if(trend_method == "OLS"){
      formula <- as.formula(paste0(dvar, "~", mvar))
      model <- lm(formula, data = A, na.action = na.omit)
    }
    
    trnd <- predict(model, B, se.fit = TRUE)
    
    if(!decreasing) {
      cdc_exc[i] <- sum(
        B[[dvar]] > trnd$fit + (conservative * sd(A[[dvar]])) &
        B[[dvar]] > (mean(A[[dvar]]) + (conservative * sd(A[[dvar]])))
      )
      cdc_p[i] <- binom.test(
        cdc_exc[i], cdc_nb[i], alternative = "greater"
      )$p.value
      cdc[i] <- if(cdc_p[i] < .05) "systematic change" else "no change"
    } else {
      cdc_exc[i] <- sum(
        B[[dvar]] < trnd$fit - (conservative * sd(A[[dvar]])) &
        B[[dvar]] < (mean(A[[dvar]]) - (conservative * sd(A[[dvar]])))
      )
      cdc_p[i] <- binom.test(
        cdc_exc[i], cdc_nb[i], alternative = "greater"
      )$p.value
      cdc[i] <- if(cdc_p[i] < .05) "systematic change" else "no change"
    }
    cdc_all <- if(length(cdc_p[cdc_p > .05]) / length(cdc_p) <= .25) {
      "systematic change"
    } else {
      "no change"
    }
  }

  out <- list(
    cdc = cdc,
    cdc_be = cdc_exc,
    cdc_b = cdc_nb,
    cdc_p = cdc_p,
    cdc_all = cdc_all,
    N = N,
    decreasing = decreasing,
    trend_method = trend_method,
    conservative = conservative,
    case_names = revise_names(data)
  )
  class(out) <- c("sc_cdc")
  attr(out, opt("phase")) <- pvar
  attr(out, opt("mt")) <- mvar
  attr(out, opt("dv")) <- dvar
  out
}