################################################################################
########################## TOPIC 1: FDR & Permutations #########################
################################################################################

load('ritalin_data.Rdata')
alpha = 0.05

# Original Data & statistics
# -------------------------------------------
yesRitalin = Data[,1:39]
noRitalin = Data[,40:78]

pValuesOriginal <- numeric(dim(Data)[1])

for (i in 1:40) {
  pValuesOriginal[i] <- t.test(yesRitalin[i,],noRitalin[i,])$p.value
}

length(which(pValuesOriginal < alpha))

# Permutations
# -------------------------------------------
nPermutations = 100

pValuesShuffled <- matrix(data = 0, nrow = dim(Data)[1], ncol = nPermutations)

for (j in 1:nPermutations) {
  permute <- sample(78)
  yesRitalinShuffled <- permute[1:39]
  noRitalinShuffled <- permute[40:78]
  
  for (i in 1:40) {
    pValuesShuffled[i,j] <- t.test(x = Data[i,yesRitalinShuffled],
                                   y = Data[i,noRitalinShuffled])$p.value
  }
}

# Computing FDR for different alphas
# -------------------------------------------
alphas <- seq(0.005,0.3,0.005)
FDR <- numeric(length(alphas))
sigOriginal<-numeric(length(alphas))

for (k in 1:length(alphas)) {
  sigOriginal[k]<- sum(pValuesOriginal < alphas[k])
  sigShuffled <- colSums(pValuesShuffled < alphas[k])
  FDR[k]<-mean(sigShuffled)/sigOriginal[k]
  FDR[k]<-max(FDR[1:k])  # avoid non-monotonic portions
}

plot(x = sigOriginal, y = FDR, type='b', xlab='# significant comparisons', ylab='FDR')


################################################################################
###################### TOPIC 2: Distances and Similarities #####################
################################################################################

load("grades.Rdata")
matplot(t(grades), type = c("b"),pch=1,xlab='Exam',ylab='Grades',lwd = 2)
legend("topright", legend = c("Student 1","Student 2","Student 3","Student 4"),col=1:4,pch=1)

# Euclidean Distance
# -----------------------------------------
BiocManager::install("ComplexHeatmap")

Euc_distance <- round(as.matrix(dist(grades,method="euclidean")),2)
Heatmap(Euc_distance, cluster_columns = F, cluster_rows = F, name = "Euclidean Distance")

# Pearson/Speaman Correlations
# -----------------------------------------
pearson_distance<-round(1-cor(t(grades),method="pearson"),2)

plot(grades[1,],grades[2,],type='p',xlim=range(grades[1,]),ylim=range(grades[2:4,]),
     pch=16,xlab="Student 1",ylab="Other students")
lines(grades[1,],grades[3,],type='p', col='blue',pch=16)
lines(grades[1,],grades[4,], type = 'p',col='red',pch=16)
legend("topright", legend = c("Student 2","Student 3","Student 4"),
       col=c("black","blue","red"),pch=16)

Heatmap(pearson_distance, cluster_columns = F, cluster_rows = F, name = "Pearson\nDistance", 
        col = c("#fff9e3","#719d92","#003f5c"))

spearman_distance<-round(1-cor(t(grades),method="spearman"),2)
Heatmap(spearman_distance, cluster_columns = F, cluster_rows = F, name = "Speaman\nDistance", 
        col = c("#fff9e3","#719d92","#003f5c"))