# 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 computing the mean of the baseline data, as well as constructing bands
# around it, according to a rule based on the standard deviation of the baseline data.
# These bands are projected into the treatment phase to explore whether the data after
# the intervention fits the pre-intervention pattern.
# The numerical output (number of points outside the bands) and the graphical representation
# may be used as a visual (although not necessarily statistical) evidence regarding the 
# presence  of change in the behavior after the intervention.
# Using this visual evidence is recommended only for baseline data not presenting 
# upward or downward trend.
#
# Suggested readings:
#
# Callahan, C. D., & Barisa, M. T. (2005). Statistical process control and rehabilitation outcome: 
# The single-subject design reconsidered. Rehabilitation Psychology, 50, 24-33. 
#
# Pfadt, A., & Wheeler, D. J. (1995). Using statistical process control to make data-based clinical decisions. 
# Journal of Applied Behavior Analysis, 28, 349-370. 
#

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

# The only part of the code that needs to be modified

# Input data
score <- c(10,8,11,6,10,6,4,5,3,4)
n_a <- 5

# Input the standard deviations' rule for creating the bands
SD_value <- 2


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


# Objects needed for the calculations
nsize <- length(score)

phaseA <- score[1:n_a]
phaseB <- score[(n_a+1):nsize]

mean_A <- rep(mean(phaseA),n_a)

# Construct bands
band <- sd(phaseA)*SD_value
upper_band <- mean(phaseA) + band
lower_band <- mean(phaseA) - band
# As vectors
upper <- rep(upper_band,nsize)
lower <- rep(lower_band,nsize)

# Counters
  
count_out_up <- 0
count_out_low <- 0
consec_out_up <- 0
consec_out_low <- 0
greatest_low <- 0
greatest_up <- 0

# Count
for (i in 1:length(phaseB))
{
  # Count below
  if (phaseB[i] < lower[n_a+i]) 
    { count_out_low <- count_out_low + 1; 
      consec_out_low <- consec_out_low + 1} else 
    { if(greatest_low < consec_out_low) greatest_low <- consec_out_low;
      consec_out_low <- 0; }
  if (greatest_low < consec_out_low) greatest_low <- consec_out_low

  
  # Count above      
  if (phaseB[i] > upper[n_a+i]) 
  { count_out_up <- count_out_up + 1; 
    consec_out_up <- consec_out_up + 1} else 
  { if(greatest_up < consec_out_up) greatest_up <- consec_out_up;
    consec_out_up <- 0}
  if (greatest_up < consec_out_up) greatest_up <- consec_out_up
}

  
# Print information
print("Number of intervention points above upper band");print(count_out_up)
print("Maximum number of consecutive intervention points above upper bands");print(greatest_up)
print("Number of intervention points below lower band");print(count_out_low)
print("Maximum number of consecutive intervention points below lower bands");print(greatest_low)

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

# Construct the plot with the SD bands

indep <- 1:nsize

# Plot limits
minimal <- min(min(score),min(lower))
maximal <- max(max(score),max(upper))

# Plot data
plot(indep,score, xlim=c(indep[1],indep[length(indep)]), ylim=c((minimal-1),(maximal+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="Standard deviation bands")

# Phase A mean & projections with SD-bands
lines(indep[1:n_a],mean_A)
lines(indep,upper,type="b")
lines(indep,lower,type="b")