#0) Create Mock Data
library(gmodels)
library(ggplot2)

#generate 'paired' data
#x1 <- round(rnorm(10, mean=50, sd=5), 1)
#x2 <- round( (x1 +rnorm(10, mean=5, sd=4)), 1)
#x3 <- round ( (x2 +rnorm(10, mean=2, sd=5)), 1)
x1 <-c(47.8, 42.0, 56.1, 52.1, 40.2, 53.1, 44.4, 57.9, 48.4, 45.3)
x2 <-c(53.6, 45.2, 61.0, 59.0, 39.2, 61.8, 48.1, 61.8, 55.2, 52.4)
x3 <-c(57.6, 46.6, 62.1, 53.8, 37.9, 63.9, 56.0, 61.1, 58.6, 52.2)

#generate datafrarme
data1 <- data.frame(value=c(x1, x2), group1=c(rep('First',10), rep('Second',10)), group2=c(rep(0,10), rep(1,10)), 
                    ID=(rep(seq(1:10),2)))
print(data1)

data2 <- data.frame(value=c(x1, x2, x3), group1=c(rep('First',10), rep('Second',10), rep('Third',10)), 
                    group2=c(rep(0,10), rep(1,10), rep(1,10)), 
                    ID=(rep(seq(1:10),3)))
print(data2)
#--------------------------------------------------------------------------------------------------------------------------


#1) Paired t-test

  #not considering paired nature
t.test(x1, x2, paired=F) #not significant
plot(data1$value~data1$group1, col="steelblue")

  #considering paired nature
t.test(x1, x2, paired=T)
plot1 <-ggplot(data=data1, aes(x=group1, y=value, group=ID)) +
  geom_line()
plot1
#--------------------------------------------------------------------------------------------------------------------------


#2) Repeated Measures ANOVA

  #not considering paired nature
model1<-aov(value~factor(group1), data=data2) #not significant
summary(model1)
plot2 <-ggplot(data=data2, aes(x=group1, y=value)) +
  geom_boxplot()
plot2

  #considering paired nature
model2<-aov(value~factor(group1)+Error(factor(ID)), data=data2)
summary(model2)

plot3 <-ggplot(data=data2, aes(x=group1, y=value, group=ID)) +
  geom_line()
plot3

  #post-hoc testing (https://www.sheffield.ac.uk/polopoly_fs/1.885219!/file/105_RepeatedANOVA.pdf)
t1 <-t.test(x1, x2, paired=T)
t2 <-t.test(x1, x3, paired=T)
t3 <-t.test(x2, x3, paired=T)
pvalues<-c(t1$p.value, t2$p.value, t3$p.value)
round(p.adjust(pvalues, 'bonferroni', 3), 4)
x   <-c(1,2,3)
y   <- c(mean(x1), mean(x2), mean(x3))
ll <- c(ci(x1)[2], ci(x2)[2], ci(x3)[2])
ul <- c(ci(x1)[3], ci(x2)[3], ci(x3)[3])
t <-c('*','','')

  #bar graph
x<-barplot(y)
br1<-barplot(y, xlab='time', ylab='measurement', ylim=c(0,70), axes=F)
xticks<-c(1,2,3)
yticks<-c(0,10,20,30,40,50,60,70)
br1 + arrows(x,ul,x,ll, length=0.05, angle=90, code=3) +box() + text(x,y+8,t) + axis(1,at=xticks) + axis(2,at=yticks)