library(agricolae)
library(multcomp)
library(emmeans)

# clustered_pearson_DS.csv - data frame for Pearson correlations
# clustered_rhoq_DS.csv - data frame for detrended correlations
# unclustered_pearson_DS.csv - data frame for Pearson correlations without eigenvector clustering
# unclustered_rhoq_DS.csv - data frame for detrended correlations without eigenvector clustering
df <- read.csv("clustered_pearson_DS.csv")
{
df$Cluster <- factor(df$Cluster)
df$Condition <- factor(df$Condition)
df$Task <- factor(df$Task)
df$Task <- relevel(df$Task,"res")
df$Part <- factor(df$Part)
levels(df$Part) <- list("rest"=c("t1_", "t2_", "t3_"), "enc"=c("enc"), "ret"=c("ret"))
df$Subject <- factor(df$Subject)
}

# Model task+condition+mode with pairwise interactions
# for 15 eigenvalues
for(i in 1:15){
{
  eval(parse(text=paste("summary(model.01 <- lm(EVal ~ (Task+Condition+Part)^2, data=df[df$Cluster==",i,",]))",sep="")))
  res=drop1(model.01, test="F")
  res.drop=which.max(res[2:dim(res)[1],6])
  # perform step-wise model reduction
  while(res[res.drop+1,6]>0.1){
    eval(parse(text=paste("summary(model.01 <- update(model.01, ~. - ",rownames(res)[res.drop+1],"))",sep="")))
    res=drop1(model.01, test="F")
    res.drop=which.max(res[2:dim(res)[1],6])
    print(rownames(res)[res.drop+1])
  }
  out <- capture.output(summary(model.01))
  cat("My title", out, file=paste("lm_evals",i,"_01.txt",sep=""), sep="\n", append=TRUE)
  summary(model.01)
}

# Do post-hoc Scheffe test
# needs library(agricolae)
{
  res=scheffe.test(model.01,"Task",group = FALSE)
  out <- capture.output(res)
  cat("My title", out, file=paste("lm_evals",i,"_Task.txt",sep=""), sep="\n", append=TRUE)
  res
  
  res=scheffe.test(model.01,"Part",group = FALSE)
  out <- capture.output(res)
  cat("My title", out, file=paste("lm_evals",i,"_Part.txt",sep=""), sep="\n", append=TRUE)
  res
}
}