# The R-code below can be used to reproduce the results from the numerical examples
# in the manuscript regarding the Douglas et al. (2015) data of patients "Samantha" and "Thomas".

# Install the packages perm, lmPerm, lmtest, nlme, SCRT, SCMA
# which are needed to perform the analyses.

install.packages("perm")
install.packages("lmPerm")
install.packages("lmtest")
install.packages("nlme")
install.packages("SCRT")
install.packages("SCMA")

# Load the packages into the workspace.

library(perm)
library(lmPerm)
library(lmtest)
library(nlme)
library(SCRT)
library(SCMA)

# Load the Samantha data into the workspace (A file selection window will open; select the Samantha.txt file).

Samantha <- read.table(file.choose(new=F), header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)

# Load the Thomas data into the workspace using the same procedure.

Thomas <- read.table(file.choose(new=F), header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)

#################################################################
# Regression analyses for single-case data of Samantha and Thomas
#################################################################

# Samantha

RegModel.1 = lm(Score~Treat, data=Samantha)
summary(RegModel.1)
confint(RegModel.1, level=0.95)

t.test(Score~Phase, alternative='two.sided', conf.level=.95, var.equal=TRUE, data=Samantha)

RegModel.1p = lmp(Score~Treat, data=Samantha)
summary(RegModel.1p)

2/choose(9,3)

permTS(Score~Phase,alternative="two.sided",exact=TRUE,method="exact.ce",data=Samantha)

RegModel.2 = lm(Score~Treat+Week, data=Samantha)
summary(RegModel.2)
confint(RegModel.2, level=0.95)

RegModel.2p = lmp(Score~Treat+Week, data=Samantha)
summary(RegModel.2p)

ArimaModel.1 <- arima(Samantha$Score, order=c(1,0,0), xreg=cbind(Samantha$Treat,Samantha$Week))
ArimaModel.1

coeftest(ArimaModel.1)

# Thomas

RegModel.1 = lm(Score~Treat, data=Thomas)
summary(RegModel.1)
confint(RegModel.1, level=0.95)

t.test(Score~Phase, alternative='two.sided', conf.level=.95, var.equal=TRUE, data=Thomas)

RegModel.1p = lmp(Score~Treat, data=Thomas)
summary(RegModel.1p)

2/choose(10,3)

permTS(Score~Phase,alternative="two.sided",exact=TRUE,method="exact.ce",data=Thomas)

RegModel.2 = lm(Score~Treat+Week, data=Thomas)
summary(RegModel.2)
confint(RegModel.2, level=0.95)

RegModel.2p = lmp(Score~Treat+Week, data=Thomas)
summary(RegModel.2p)

ArimaModel.1 <- arima(Thomas$Score, order=c(1,0,0), xreg=cbind(Thomas$Treat,Thomas$Week))
ArimaModel.1

coeftest(ArimaModel.1)

#################################################
# Multilevel analysis of Samantha and Thomas data
#################################################

Total <- read.table(file.choose(new=F), header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
# Select the file: "Combined_data.txt"

Model.0 <- lme(Score~Treat, random=~1|Case, data=Total, control=list(opt="optim"))
summary(Model.0)
VarCorr(Model.0)
intervals(Model.0)
varests <- as.numeric(VarCorr(Model.0)[1:2])
ICC <- varests[1]/sum(varests)
ICC
sqrt(diag(Model.0$apVar))

Model.1 <- lme(Score~Treat, random=~Treat|Case, data=Total, control=list(opt="optim"))
# optim problem, convergence error code = 1
summary(Model.1)
VarCorr(Model.1)
intervals(Model.1)
varests <- as.numeric(VarCorr(Model.1)[1:2])
ICC <- varests[1]/sum(varests)
ICC

Model.1b <- lme(Score~Treat, random = list(Case = pdDiag(~Treat)), data=Total, control=list(opt="optim"))
summary(Model.1b)
VarCorr(Model.1b)
intervals(Model.1b)
varests <- as.numeric(VarCorr(Model.1b)[1:2])
ICC <- varests[1]/sum(varests)
ICC
sqrt(diag(Model.1b$apVar))

Model.2 <- lme(Score~Treat+Time, random=~Time+Treat|Case, data=Total, control=list(opt="optim"))
#optim problem, convergence error code = 1
summary(Model.2)
VarCorr(Model.2)
intervals(Model.2)
varests <- as.numeric(VarCorr(Model.2)[1:2])
ICC <- varests[1]/sum(varests)
ICC

Model.2b <- lme(Score~Treat+Time, random = list(Case = pdDiag(~Treat+Time)), data=Total, control=list(opt="optim"))
summary(Model.2b)
VarCorr(Model.2b)
intervals(Model.2b)
varests <- as.numeric(VarCorr(Model.2b)[1:2])
ICC <- varests[1]/sum(varests)
ICC
sqrt(diag(Model.2b$apVar))

Model.3b <- update(Model.2b, correlation=corAR1())
summary(Model.3b)
VarCorr(Model.3b)
anova(Model.2b,Model.3b)
sqrt(diag(Model.2b$apVar))

############################################################################
# Randomization tests for the Samantha and Thomas data and P-value combining
############################################################################

pvalue.systematic(design = "AB", statistic = "|A-B|", limit = 3, data = Samantha)
pvalue.systematic(design = "AB", statistic = "|A-B|", limit = 3, data = Thomas)

Pvalues <- read.table(file.choose(new=F), header=TRUE, sep="", na.strings="NA", dec=".", strip.white=TRUE)
# Select file: "Pvalues.txt"

combine(method = "+", pvalues = Pvalues)

(.45)*(.45)/2

distribution.systematic(design = "AB", statistic = "|A-B|", limit = 3, data = Samantha)
observed(design = "AB", statistic = "|A-B|", data = Samantha)

################################################################################
# Randomization test wrapper for multilevel analysis of Samantha and Thomas data
################################################################################

# N Samantha=9
# N Thomas=10

# Step 1: Construct random startpoint assignment matrix for the Samantha and Thomas data.

quantity(MT=9, design="AB", limit=3)
quantity(MT=10, design="AB", limit=3)

ass_S <- assignments(design="AB", MT=9, limit=3) # all possible assignments for Samantha
ass_T <- assignments(design="AB", MT=10, limit=3) # all possible assignments for Thomas

ass_combined <- matrix(nrow=19, ncol=20)

c <- 1
for (i in 1:4)
{
  for(j in 1:5)
  {
   ass_combined[,c] <- c(ass_S[i,], ass_T[j,])
   c <- c+1
  }
}

ass_combined # combined random startpoint assignment matrix for the Samantha and Thomas data (1 column = 1 assignment) 

t_values <- numeric()
p_values <- numeric()

# Step 2: Perform multilevel analysis for each design of the random startpoint assignment matrix
# and record the resulting t-value and p-value

for (i in 1:20)
{
  Total_r <- Total
  Total_r[,5] <- ass_combined[,i]
  
  Model.2b <- lme(Score~Treat+Time, random = list(Case = pdDiag(~Treat+Time)), data=Total_r, control=list(opt="optim"))
  t_values[i] <- summary(Model.2b)$tTable[2,"t-value"] # tvalue
  p_values[i] <- summary(Model.2b)$tTable[2,"p-value"] # pvalue
}

# Step 3: Give a name to each possible design and construct a dataframe that
# displays the t- and p-values corresponding to each design


designs <- c("(4;4)","(4;5)","(4;6)","(4;7)","(4;8)",
  "(5;4)","(5;5)","(5;6)","(5;7)","(5;8)",
  "(6;4)","(6;5)","(6;6)","(6;7)","(6;8)",
  "(7;4)","(7;5)","(7;6)","(7;7)","(7;8)")

results <- data.frame(designs,t_values, p_values)

# Order the dataframe according to decreasing t-value

results[order(t_values, decreasing=T),]