# R SYNTAX COPYRIGHT (C) 2013 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 phase difference" procedure
# described in Manolov and Solanas (2013):
# 
# Manolov, R., & Solanas, A. (2013). A comparison of mean phase difference 
#   and generalized least squares for analyzing single-case data. 
#   Journal of School Psychology, 51, 201-215.
#
# The R script also offers a graphical representation of the actual and projected data.


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

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


# Example data set: baseline measurements change the values within ()
baseline <- c(1,2,3)
# Example data set: intervention phase measurements change the values within ()
treatment <- c(5,6,7)

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

# THE FOLLOWING CODE NEEDS NOT BE CHANGED

# Obtain phase length
n_a <- length(baseline)
n_b <- length(treatment)

# Estimate baseline trend
base_diff <- rep(0,(n_a-1))
for (i in 1:(n_a-1))
    base_diff[i] <- baseline[i+1] - baseline[i]
trendA <- mean(base_diff)

# Project baseline trend to the treatment phase
treatment_pred <- rep(0,n_b)
for (i in 1:n_b) 
    treatment_pred[i] <- baseline[1] + trendA*(n_a+i-1)

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

# Represent graphically the actual data and the projected treatment data

info <- c(baseline,treatment)
info_pred <- c(baseline,treatment_pred)
minimo <- min(info,info_pred)
maximo <- max(info,info_pred)
time <- c(1:length(info))
plot(time,info, xlim=c(1,length(info)), ylim=c((minimo-1),(maximo+1)), xlab="Measurement time", ylab="Behavior 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):length(info)],info[(n_a+1):length(info)])
axis(side=1, at=seq(0,length(info),1),labels=TRUE, font=2)
points(time, info, pch=24, bg="black")
points (time, info_pred, pch=19)
title (main="Predicted (circle) vs. actual (triangle) data")

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

# PERFORM THE CALCULATIONS AND PRINT THE RESULT

# Compare the predicted and the actually obtained treatment data and obtain the average difference
mpd <- mean(treatment-treatment_pred)
print ("The mean phase difference is"); print(mpd)