# BASELINE CORRECTED TAU FOR SINGLE-CASE RESEARCH (2017, JUNE)

# Citation: Tarlow, K. R. (2017, June). Baseline Corrected Tau
#           for single-case research (R code). Retrieved from 
#           http://ktarlow.com/stats/

# R SYNTAX COPYRIGHT (C) 2017 KEVIN R. TARLOW
# url: http://www.ktarlow.com/stats
# email: krtarlow@gmail.com

# adapted from: Tarlow, K. R. (2017). An improved rank correlation
#               effect size statistic for single-case designs: 
#               Baseline Corrected Tau. Behavior Modification, 41(4), 
#               427-467.

# This work is licensed under the Creative Commons 
# Attribution-NonCommercial 3.0 Unported License. 
# To view a copy of this license, visit 
# http://creativecommons.org/licenses/by-nc/3.0/deed.en_US
# 
# You are free to copy, distribute, transmit, and adapt
# the work under the following conditions:
#
# Attribution - You must attribute the work in the manner 
# specified by the author, KEVIN R. TARLOW (but not in any way 
# that suggests that the author endorses you or your use 
# of the work).
#
# Noncommercial  You may not use this work for commercial 
# purposes.

theil <- function(x,y) {

# returns theil-sen slope and intercept estimates;
# x and y are two equal length vectors (x & y coords)

	n <- length(x)
	slopes <- as.numeric()
	ints <- as.numeric()
	for (i in 1:(n-1)) {
		for(j in (i+1):n) {
			slopes <- c(slopes, ((y[j] - y[i]) / (x[j] - x[i])))
		}
	}
	b <- median(slopes)
	for (i in 1:n) {
		ints <- c(ints, (y[i] - (b*x[i])))
	}
	results <- list(slope=b, int=median(ints))
	return(results)
}

bctau <- function(a,b) {

# The bctau() function accepts two arguments, a and b, which
# are vectors for each phase in an AB single-case design

	library(Kendall)

	n <- length(a) + length(b)

	ta <- 1:(length(a))
	tb <- (length(a)+1):(length(a)+length(b))

	# if baseline trend is not statistically significant, 
	# return tau result (no trend correction)

		if(Kendall(a,ta)$sl > .05) {
			results <- Kendall(c(a,b),c(rep(0,length(a)),rep(1,length(b))))
			tau <- as.numeric(results$tau)
			p <- as.numeric(results$sl)
			se <- sqrt((2/n) * (1-(tau^2)))
			return(list(tau=tau, p=p, se=se, corrected=FALSE))
		}

	# if baseline trend is statistically significant,
	# get Theil-Sen residuals

		theilsen <- theil(ta, a)
		slope <- theilsen$slope
		intercept <- theilsen$int

		correcteda <- as.numeric()
		correctedb <- as.numeric()

		for (i in 1:length(a)) {
			correcteda[i] <- a[i] - (slope*i + intercept)
		}
		for (i in 1:length(b)) {
			correctedb[i] <- b[i] - (slope*(i+length(a)) + intercept)
		}

		results <- Kendall(c(correcteda,correctedb),c(rep(0,length(a)),rep(1,length(b))))
		tau <- as.numeric(results$tau)
		p <- as.numeric(results$sl)
		se <- sqrt((2/n) * (1-(tau^2)))
		return(list(tau=tau, p=p, se=se, corrected=TRUE, int=intercept, slope=slope, correcteda=correcteda, correctedb=correctedb))

}