# R SYNTAX COPYRIGHT (C) 2010 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 "Slope and level change" procedure
# described in Solanas, Manolov, and Onghena (2010):
# 
# Solanas, A., Manolov, R., & Onghena, P. (2010). Estimating slope and level change 
#   in N=1 designs. Behavior Modification, 34, 195-218.
#
# The R script also offers a graphical representation of the actual and detrended data.


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

# MODIFY THE EXAMPLE AB-DATA SET ACCORDING TO YOUR DATA

# Example data set AB data: change the values within ()
info <- c(10,9,9,8,7,8,8,7,7,6,4,3,2,1,1,2,1,1,0,0)
# Baseline phase length: change the value according to your data
n_a <- 10

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

# THE FOLLOWING CODE NEEDS NOT BE CHANGED

# Initial data manipulations
slength <- length(info)
n_b <- slength-n_a
phaseA <- info[1:n_a]
phaseB <- info[(n_a+1):slength]

# Estimate phase A trend
phaseAdiff <- c(1:(n_a-1))
for (iter in 1:(n_a-1))
  phaseAdiff[iter] <- phaseA[iter+1] - phaseA[iter]
trendA <- mean(phaseAdiff)

# Remove phase A trend from the whole data series
phaseAdet <- c(1:n_a)
for (timeA in 1:n_a)
  phaseAdet[timeA] <- phaseA[timeA] - trendA * timeA

phaseBdet <- c(1:n_b)
for (timeB in 1:n_b)
  phaseBdet[timeB] <- phaseB[timeB] - trendA * (timeB+n_a)

# Compute the slope change estimate
phaseBdiff <- c(1:(n_b-1))
for (iter in 1:(n_b-1))
  phaseBdiff[iter] <- phaseBdet[iter+1] - phaseBdet[iter]

trendB <- mean(phaseBdiff)


# Compute the level change estimate
phaseBddet <- c(1:n_b)
for (timeB in 1:n_b)
  phaseBddet[timeB] <- phaseBdet[timeB] - trendB * (timeB-1)

level <- mean(phaseBddet) - mean(phaseAdet)



######################################################################
# Represent graphically actual and detrended data

time <- c(1:slength)
par(mfrow=c(2,1))
plot(time,info, xlim=c(1,slength), ylim=c((min(info)-1),(max(info)+1)), xlab="Measurement time", ylab="Variable of interest", font.lab=2)
abline(v=(n_a+0.5))
lines(time[1:n_a],info[1:n_a])
lines(time[(n_a+1):slength],info[(n_a+1):slength])
axis(side=1, at=seq(0,slength,1),labels=TRUE, font=2)
axis(side=2, at=seq((min(info)-1),(max(info)+1),2),labels=TRUE, font=2)
points(time, info, pch=24, bg="black")
title (main="Original data")
transf <- c(phaseAdet,phaseBdet)
plot(time,transf, xlim=c(1,slength), ylim=c((min(transf)-1),(max(transf)+1)), xlab="Measurement time", ylab="Variable of interest", font.lab=2)
abline(v=(n_a+0.5))
lines(time[1:n_a],transf[1:n_a])
lines(time[(n_a+1):slength],transf[(n_a+1):slength])
axis(side=1, at=seq(0,slength,1),labels=TRUE, font=2)
axis(side=2, at=seq((min(transf)-1),(max(transf)+1),2),labels=TRUE, font=2)
points(time, transf, pch=24, bg="black")
title (main="Detrended data")

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

# PRINT THE RESULT

print ("Slope change estimate = "); print(trendB)
print ("Net level change estimate = "); print(level)