# R SYNTAX COPYRIGHT (C) 2009 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 "Percentage of nonoverlapping data" index
# described in Manolov and Solanas (2009):
# 
# Manolov, R., & Solanas, A. (2009). Percentage of nonoverlapping corrected data. 
#   Behavior Research Methods, 41, 1262-1271.
#
# 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: baseline measurements change the values within ()
phaseA <- c(9,8,8,7,6,7,7,6,6,5)
# Example data set: intervention phase measurements change the values within ()
phaseB <- c(5,6,4,3,3,6,2,2,2,1)

# Specify whether the intervention should increase or reduce the behavior
aim <- "reduce" # Alternatively aim <- "increase"

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

# THE FOLLOWING CODE NEEDS NOT BE CHANGED

n_a <- length(phaseA)
n_b <- length(phaseB)

# Data correction: phase A
phaseAdiff <- c(1:(n_a-1))
for (iter1 in 1:(n_a-1))
phaseAdiff[iter1] <- phaseA[iter1+1] - phaseA[iter1]
phaseAcorr <- c(1:n_a)
for (iter2 in 1:n_a)
phaseAcorr[iter2] <- phaseA[iter2] - mean(phaseAdiff)*iter2

# Data correction: phase B
phaseBcorr <- c(1:n_b)
for (iter3 in 1:n_b)
phaseBcorr[iter3] <- phaseB[iter3] - mean(phaseAdiff)*(iter3+n_a)


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

slength <- n_a + n_b
info <- c(phaseA,phaseB)
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(phaseAcorr,phaseBcorr)
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")

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

# PERFORM THE CALCULATIONS AND PRINT THE RESULT

# PND on corrected data: Aim to increase
if (aim == "increase")
{ 
  countcorr <- 0
  for (iter4 in 1:n_b)
    if (phaseBcorr[iter4] > max(phaseAcorr)) countcorr <- countcorr+1

  pndcorr <- (countcorr/n_b)*100
  print ("The percent of nonoverlapping corrected data is"); print(pndcorr)
}


# PND on corrected data: Aim to reduce
if (aim == "reduce")
{ 
  countcorr <- 0
  for (iter4 in 1:n_b)
    if (phaseBcorr[iter4] < min(phaseAcorr)) countcorr <- countcorr+1

  pndcorr <- (countcorr/n_b)*100
  print ("The percent of nonoverlapping corrected data is"); print(pndcorr)
}