# R SYNTAX COPYRIGHT (C) 2014 Rumen Manolov
#
# 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, Rumen Manolov (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.
#

# This R script allows estimating baseline trend via the split-middle method
# (Miller, 1985), projecting the trend into the treatment phase and constructing
# two types of envelope around this projection: on the basis of the baseline phase
# median (Gast & Spriggs, 2010) and on the basis of the interquartile range (Tukey, 1977).
#
#Gast, D. L., & Spriggs, A. D. (2010). Visual analysis of graphic data. In D. L. Gast (Ed.), 
#Single subject research methodology in behavioral sciences (pp. 199-233). London, UK: Routledge.
#
# Miller, M. J. (1985). Analyzing client change graphically. Journal of Counseling and Development, 63, 491-494. 
#
# Tukey, J. W. (1977). Exploratory data analysis. London, UK: Addison-Wesley.
#
# The code accompanies the following study:
# 
# Manolov, R., Sierra, V., Solanas, A., & Botella, J. (2014). Assessing functional relations in single-case designs: 
# Quantitative proposals in the context of the evidence-based movement. Unpublished manuscript.

# The R script also offers a graphical representation of the baselie split-middle trend and the upper and lower limits of its projection
# It also offers a quantification of the proportion of treatment phase scores falling within the two types of envelope
# constructed around the projected trend line.

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

# Input data
score <- c(1,3,5,5,10,9,11,14)
n_a <- 4
# Input the percentage of the Median to use for constructing the envelope
md_percentage <- 20
# Input the interquartile range to use for constructing the envelope
IQR_value <- 1.5
# Choose figures display
display <- "vertical" # Alternatively "horizontal"

#####################################################
# This part of the code needs not be changed: only copy-paste it in the R console

# Objects needed for the calculations
md_addsub <- md_percentage/200
nsize <- length(score)
indep <- 1:nsize

# Construct phase vectors
a1 <- score[1:n_a]
b1 <- score[(n_a+1):nsize]
 
 
##############################
# Split-middle method

# Divide the phase into two parts of equal size: even number of measm
  if (length(a1)%%2==0) 
  { 
    part1 <- 1:(length(a1)/2)
    part1 <- a1[1:(length(a1)/2)]
    part2 <- ((length(a1)/2)+1):length(a1)
    part2 <- a1[((length(a1)/2)+1):length(a1)]
  } 


  # Divide the phase into two parts of equal size: odd number of measm
  if (length(a1)%%2==1) 
  { 
    part1 <- 1:((length(a1)-1)/2)
    part1 <- a1[1:((length(a1)-1)/2)]
    part2 <- (((length(a1)+1)/2)+1):length(a1)
    part2 <- a1[(((length(a1)+1)/2)+1):length(a1)]
  } 

  # Obtain the median in each section
  median1 <- median(part1)
  median2 <- median(part2)

  # Obtain the midpoint in each section
  midp1 <- (length(part1)+1)/2
  if (length(a1)%%2==0) midp2 <- length(part1) + (length(part2)+1)/2
  if (length(a1)%%2==1) midp2 <- length(part1) + 1 + (length(part2)+1)/2

  # Obtain the values of the split-middle trend: "regla de 3"
  slope <- (median2-median1)/(midp2-midp1)
  sm.trend <- 1:length(a1)


  # Whole number function
  is.wholenumber <-
    function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

  if (!is.wholenumber(midp1)) 
  { 
    sm.trend[midp1-0.5] <- median1 - (slope/2)
    for (i in 1:length(part1))
      if ((midp1-0.5 - i) >= 1) sm.trend[(midp1-0.5 - i)] <- sm.trend[(midp1-0.5-i + 1)] - slope
    for (i in (midp1+0.5):length(sm.trend))
      sm.trend[i] <- sm.trend[i-1] + slope 
  }

  if (is.wholenumber(midp1)) 
  {
    sm.trend[midp1] <- median1
    for (i in 1:length(part1))
      if ((midp1 - i) >= 1) sm.trend[(midp1 - i)] <- sm.trend[(midp1-i + 1)] - slope
    for (i in (midp1+1):length(sm.trend))
      sm.trend[i] <- sm.trend[i-1] + slope
  }

# Consutructing envelopes
 
  # Construct envelope 20% median (Gast & Spriggs, 2010)
  envelope <- median(a1)*md_addsub
  upper_env <- 1:length(a1) 
  upper_env <- sm.trend + envelope
  lower_env <- 1:length(a1) 
  lower_env <- sm.trend - envelope

  # Construct envelope 1.5 IQR (add to median)
  tukey.eda <- IQR_value*IQR(a1)
  upper_IQR <- 1:length(a1) 
  upper_IQR <- sm.trend + tukey.eda
  lower_IQR <- 1:length(a1) 
  lower_IQR <- sm.trend - tukey.eda

  # Project split middle trend and envelope
  projected <- 1:length(b1)
  proj_env_U <- 1:length(b1)
  proj_env_L <- 1:length(b1)
  proj_IQR_U <- 1:length(b1)
  proj_IQR_L <- 1:length(b1)

  projected[1] <- sm.trend[length(sm.trend)]+slope
  proj_env_U[1] <- upper_env[length(upper_env)]+slope
  proj_env_L[1] <- lower_env[length(lower_env)]+slope
  proj_IQR_U[1] <- upper_IQR[length(upper_IQR)]+slope
  proj_IQR_L[1] <- lower_IQR[length(lower_IQR)]+slope

  for (i in 2:length(b1))
  {
    projected[i] <- projected[i-1]+slope
    proj_env_U[i] <- proj_env_U[i-1]+slope
    proj_env_L[i] <- proj_env_L[i-1]+slope
    proj_IQR_U[i] <- proj_IQR_U[i-1]+slope
    proj_IQR_L[i] <- proj_IQR_L[i-1]+slope
  }

  # Count the number of values within the envelope
  in_env <-0
  in_IQR <- 0

  for (i in 1:length(b1))
  {
    if ( (b1[i] >= proj_env_L[i]) && (b1[i] <= proj_env_U[i]) ) in_env <- in_env + 1
    if ( (b1[i] >= proj_IQR_L[i]) && (b1[i] <= proj_IQR_U[i]) ) in_IQR <- in_IQR + 1
  }  

  prop_env <- in_env/length(b1)
  prop_IQR <- in_IQR/length(b1)

  
# Print information
print("Proportion of phase B data into envelope");print(prop_env)
print("Proportion of phase B data into IQR limits");print(prop_IQR)


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

# Construct the plot with the median-based envelope

if (display=="vertical") {rows <- 2; cols <- 1} else  {rows <- 1; cols <- 2}

minimal1 <- min(min(score),min(proj_env_L))
maximal1 <- max(max(score),max(proj_env_U))

par(mfrow=c(rows,cols))

# Plot data
plot(indep,score, xlim=c(indep[1],indep[length(indep)]), ylim=c((minimal1-1),(maximal1+1)), xlab="Measurement time", ylab="Score", font.lab=2)
lines(indep[1:n_a],score[1:n_a])
lines(indep[(n_a+1):nsize],score[(n_a+1):nsize])
abline (v=(n_a+0.5))
points(indep, score, pch=24, bg="black")
title(main="Median-based envelope around projected split-middle trend")

# Split-middle trend & projections with limits
lines(indep[1:n_a],sm.trend[1:n_a])
lines(indep[(n_a+1):nsize],proj_env_U,type="b")
lines(indep[(n_a+1):nsize],proj_env_L,type="b")


# Construct the plot with the IQR-based envelope

minimal2 <- min(min(score),min(proj_IQR_L))
maximal2 <- max(max(score),max(proj_IQR_U))

# Plot data
plot(indep,score, xlim=c(indep[1],indep[length(indep)]),ylim=c((minimal2-1),(maximal2+1)), xlab="Measurement time", ylab="Score", font.lab=2)
lines(indep[1:n_a],score[1:n_a])
lines(indep[(n_a+1):nsize],score[(n_a+1):nsize])
abline (v=(n_a+0.5))
points(indep, score, pch=24, bg="black")
title(main="IQR-based envelope around projected split-middle trend")

# Split-middle trend & projections with limits
lines(indep[1:n_a],sm.trend[1:n_a])
lines(indep[(n_a+1):nsize],proj_IQR_U,type="b")
lines(indep[(n_a+1):nsize],proj_IQR_L,type="b")