"
     _    _ _____  _____ _____ 
    | |  | /  __ \/  ___|_   _|
    | |  | | /  \/\ `--.  | |  
    | |/\| | |     `--. \ | |  
    \  /\  / \__/\/\__/ / | |  
     \/  \/ \____/\____/  \_/  
  
      The Wisconsin Project
"
#----------------------------------------------------------------------------------#

      # Complete R syntax
      # by Philip Lindner (philip.sebastian.lindner@gmail.com)
      # Stockholm University, October 2015

#----------------------------------------------------------------------------------#

  # Loading data
  old.par <- par()
  wcst <- read.csv2("Projects_current/WCST/data.csv")

  # Excluding bad cases (explained in manuscript) and naming groups 
  good.wcst <- subset(wcst, wcst$exclude == 0)
  good.wcst$study <- factor(good.wcst$study,levels=c(1,2,3),labels=c("SAD","TIN","MDD"))
  
  # Creating variable raw score difference between pre- and post-treatment symptom measures, and percentage PE of TE
  good.wcst$raw_pom_diff <- good.wcst$raw_pom_pre - good.wcst$raw_pom_post
  good.wcst$percPE <- good.wcst$PE / good.wcst$TE
  
  # Partialling out age out of WCST metrics (akin to using age-adjusted norm data)
  lm.PE.age <- lm(good.wcst$PE ~ good.wcst$age)
  lm.CC.age <- lm(good.wcst$CC ~ good.wcst$age)
  lm.TC.age <- lm(good.wcst$TC ~ good.wcst$age)
  lm.T1C.age <- lm(good.wcst$T1C ~ good.wcst$age)
  lm.SEC.age <- lm(good.wcst$SEC ~ good.wcst$age)
  lm.percPE.age <- lm(good.wcst$percPE ~ good.wcst$age)
  lm.TE.age <- lm(good.wcst$TE ~ good.wcst$age)
  good.wcst$adjPE <- residuals(lm.PE.age)
  good.wcst$adjCC <- residuals(lm.CC.age)
  good.wcst$adjTC<- residuals(lm.TC.age)
  good.wcst$adjT1C <- residuals(lm.T1C.age)
  good.wcst$adjSEC <- residuals(lm.SEC.age)
  good.wcst$adjpercPE <- residuals(lm.percPE.age)
  good.wcst$adjTE <- residuals(lm.TE.age)
  
  View(good.wcst)
  
  # Subsetting
  good.treatment.wcst <- subset(good.wcst, any_beh == 1)
  
  good.sad.wcst <- subset(good.wcst,study == "SAD")
  good.sad.treatment.wcst <- subset(good.sad.wcst, any_beh == 1)
    
  good.tin.wcst <- subset(good.wcst,study == "TIN")
  good.tin.treatment.wcst <- subset(good.tin.wcst, any_beh == 1)
  
  good.mdd.wcst <- subset(good.wcst,study == "MDD")
  good.mdd.treatment.wcst <- subset(good.mdd.wcst, any_beh == 1)
  
  
### Sample descriptions mentioned in manuscript
  
  library(psych)
  
  describe(good.wcst$SEC)

  table(good.wcst$study)
  table(good.wcst$any_beh,good.wcst$study)
  
  # Table 1
  
  describe(good.wcst$age)
  describeBy(good.wcst$age,good.wcst$study,range=F,skew=F)
  summary(aov(age ~ study, data=good.wcst))
  pairwise.t.test(good.wcst$age,good.wcst$study,p.adj="bonf")
  
  t1.temp <- table(good.wcst$sex, good.wcst$study)
  100*round(prop.table(t1.temp, 2),3)
  fisher.test(good.wcst$sex, good.wcst$study)["p.value"]
      # No sex-difference in any WCST metric (not in Table 1 but in Methods section)
      summary(aov(PE ~ sex, good.wcst))
      summary(aov(TE ~ sex, good.wcst))
      summary(aov(CC ~ sex, good.wcst))
      summary(aov(TC ~ sex, good.wcst))
      summary(aov(T1C ~ sex, good.wcst))
      summary(aov(SEC ~ sex, good.wcst))
  
  # Raw scores
  
  summary(aov(PE ~ study, good.wcst))
  describeBy(good.wcst$PE, good.wcst$study,range=F,skew=F)
  
  summary(aov(TE ~ study, good.wcst))
  describeBy(good.wcst$TE, good.wcst$study,range=F,skew=F)
  
  summary(aov(CC ~ study, good.wcst))
  describeBy(good.wcst$CC, good.wcst$study,range=F,skew=F)
  
  summary(aov(TC ~ study, good.wcst))
  describeBy(good.wcst$TC, good.wcst$study,range=F,skew=F)
  
  summary(aov(T1C ~ study, good.wcst))
  describeBy(good.wcst$T1C, good.wcst$study,range=F,skew=F)
  
  summary(aov(SEC ~ study, good.wcst))
  describeBy(good.wcst$SEC, good.wcst$study,range=F,skew=F)
  pairwise.t.test(good.treatment.wcst$SEC, good.treatment.wcst$study,p.adj="bonf")
  
  # Age-adjusted
  
  summary(aov(adjPE ~ study, good.wcst))
  describeBy(good.wcst$adjPE, good.wcst$study,range=F,skew=F)
  
  summary(aov(adjTE ~ study, good.wcst))
  describeBy(good.wcst$adjTE, good.wcst$study,range=F,skew=F)
  
  summary(aov(adjCC ~ study, good.wcst))
  describeBy(good.wcst$adjCC, good.wcst$study,range=F,skew=F)
  
  summary(aov(adjTC ~ study, good.wcst))
  describeBy(good.wcst$adjTC, good.wcst$study,range=F,skew=F)
  
  summary(aov(adjT1C ~ study, good.wcst))
  describeBy(good.wcst$adjT1C, good.wcst$study,range=F,skew=F)
  
  summary(aov(adjSEC ~ study, good.wcst))
  describeBy(good.wcst$adjSEC, good.wcst$study,range=F,skew=F)
  
  summary(aov(PSR ~ study, good.treatment.wcst))
  pairwise.t.test(good.treatment.wcst$PSR,good.treatment.wcst$study,p.adj="bonf")
  

### Factor analysis
  
  
  forfactoranalysis <- c("CC","TC","PE","TE")
  good.wcst.fa <- good.wcst[forfactoranalysis]
  principal(good.wcst.fa, nfactors=2, rotate="none")
  
    
### Results 3.1. Linear associations

  # Pre-treatment scores
  
  lm.pre.sad <- (lm(raw_pom_pre ~ adjPE,good.sad.wcst))
  summary(lm.pre.sad)
  round(confint(lm.pre.sad)[c(2,4)],3)
  lm.pre.sad.intercept <- coef(lm.pre.sad)[1]
  lm.pre.sad.pe <- coef(lm.pre.sad)[2]
  
  lm.pre.tin <- (lm(raw_pom_pre ~ adjPE, good.tin.wcst))
  summary(lm.pre.tin)
  round(confint(lm.pre.tin)[c(2,4)],3)
  lm.pre.tin.intercept <- coef(lm.pre.tin)[1]
  lm.pre.tin.pe <- coef(lm.pre.tin)[2]
  
  lm.pre.mdd <- (lm(raw_pom_pre ~ adjPE,good.mdd.wcst))
  summary(lm.pre.mdd)
  round(confint(lm.pre.mdd)[c(2,4)],3)
  lm.pre.mdd.intercept <- coef(lm.pre.mdd)[1]
  lm.pre.mdd.pe <- coef(lm.pre.mdd)[2]
  
  # Raw outcome measures
  
  lm.pom.sad <- (lm(raw_pom_diff ~ adjPE,good.sad.treatment.wcst))
  summary(lm.pom.sad)
  round(confint(lm.pom.sad)[c(2,4)],3)
  lm.pom.sad.intercept <- coef(lm.pom.sad)[1]
  lm.pom.sad.pe <- coef(lm.pom.sad)[2]
  
  lm.pom.tin <- (lm(raw_pom_diff ~ adjPE, good.tin.treatment.wcst))
  summary(lm.pom.tin)
  round(confint(lm.pom.tin)[c(2,4)],3)
  lm.pom.tin.intercept <- coef(lm.pom.tin)[1]
  lm.pom.tin.pe <- coef(lm.pom.tin)[2]
  
  lm.pom.mdd <- (lm(raw_pom_diff ~ adjPE,good.mdd.treatment.wcst))
  summary(lm.pom.mdd)
  round(confint(lm.pom.mdd)[c(2,4)],3)
  lm.pom.mdd.intercept <- coef(lm.pom.mdd)[1]
  lm.pom.mdd.pe <- coef(lm.pom.mdd)[2]
  
  # PSR

  lm.psr.all <- (lm(PSR ~ adjPE,good.treatment.wcst))
  summary(lm.psr.all)
  round(confint(lm.psr.all)[c(2,4)],3)
  lm.psr.all.intercept <- coef(lm.psr.all)[1]
  lm.psr.all.pe <- coef(lm.psr.all)[2]

  lm.psr.sad <- (lm(PSR ~ adjPE,good.sad.treatment.wcst))
  summary(lm.psr.sad)
  summary(lm(PSR ~ adjPE,good.sad.treatment.wcst[which(good.sad.treatment.wcst$PSR>-50),] )) #Re-running excluding PSR-outliers
  round(confint(lm.psr.sad)[c(2,4)],3)
  lm.psr.sad.intercept <- coef(lm.psr.sad)[1]
  lm.psr.sad.pe <- coef(lm.psr.sad)[2]
  
  lm.psr.tin <- (lm(PSR ~ adjPE, good.tin.treatment.wcst))
  summary(lm.psr.tin)
  round(confint(lm.psr.tin)[c(2,4)],3)
  lm.psr.tin.intercept <- coef(lm.psr.tin)[1]
  lm.psr.tin.pe <- coef(lm.psr.tin)[2]
  
  lm.psr.mdd <- (lm(PSR ~ adjPE,good.mdd.treatment.wcst))
  summary(lm.psr.mdd)
  round(confint(lm.psr.mdd)[c(2,4)],3)
  lm.psr.mdd.intercept <- coef(lm.psr.mdd)[1]
  lm.psr.mdd.pe <- coef(lm.psr.mdd)[2]
  
  #### Percentage PE 
    
  # Pre-treatment scores
  summary(lm(raw_pom_pre ~ adjpercPE,good.sad.wcst))
  summary(lm(raw_pom_pre ~ adjpercPE, good.tin.wcst))
  summary(lm(raw_pom_pre ~ adjpercPE,good.mdd.wcst))
  
    
  # Raw outcome measures
  summary(lm(raw_pom_diff ~ adjpercPE,good.sad.treatment.wcst))
  summary(lm(raw_pom_diff ~ adjpercPE, good.tin.treatment.wcst))
  summary(lm(raw_pom_diff ~ adjpercPE,good.mdd.treatment.wcst))
  
  # PSR
  summary(lm(PSR ~ adjpercPE,good.treatment.wcst))
  summary(lm(PSR ~ adjpercPE,good.sad.treatment.wcst))
  summary(lm(PSR ~ adjpercPE, good.tin.treatment.wcst))
  summary(lm(PSR ~ adjpercPE,good.mdd.treatment.wcst))

  

  
#### Figure 2 -- Scatterplots (export as 14*14 PDF)
  
  # Loading graphics packages
  
  library(ggplot2)
  library(gridExtra)
  
  plot.all <- ggplot(good.treatment.wcst, aes(x=PE,y=PSR,color=study)) + expand_limits(x=c(0,40)) + theme(legend.position="none") + ggtitle("Full sample and percentage symptom reduction (PSR)") + geom_point(size=4) + scale_color_manual(values=c("coral","chartreuse3","cornflowerblue")) + ylab("Percentage symptom reduction") + xlab("Perseverative errors (raw count)") 
  
  plot.sad <- ggplot(good.sad.treatment.wcst,aes(x=PE,y=raw_pom_diff)) + expand_limits(x=c(0,40)) + geom_point(size=4,color="coral") + ggtitle("Social anxiety disorder sample and LSAS-SR pre-post score difference") + ylab("LSAS-SR pre-post score difference") + xlab("Perseverative errors (raw count)")
  
  plot.tin <- ggplot(good.tin.treatment.wcst,aes(x=PE,y=raw_pom_diff)) + expand_limits(x=c(0,40)) + geom_point(size=4,color="chartreuse3") + ggtitle("Tinnitus sample and THI pre-post score difference")  + ylab("THI pre-post score difference") + xlab("Perseverative errors (raw count)")
  
  plot.mdd <- ggplot(good.mdd.treatment.wcst,aes(x=PE,y=raw_pom_diff)) + expand_limits(x=c(0,40)) + geom_point(size=4,color="cornflowerblue") + ggtitle("Depression sample and BDI-II pre-post score difference") + ylab("BDI-II pre-post score difference") + xlab("Perseverative errors (raw count)")
    
  grid.arrange(plot.all, plot.sad,plot.tin,plot.mdd,nrow=2)
  
      # Regression lines used in initial PeerJ submission but since omitted
      + geom_abline(intercept=lm.psr.all.intercept, slope=lm.psr.all.pe) + geom_abline(intercept=lm.psr.sad.intercept, slope=lm.psr.sad.pe, color="coral") + geom_abline(intercept=lm.psr.tin.intercept, slope=lm.psr.tin.pe, color="chartreuse3") + geom_abline(intercept=lm.psr.mdd.intercept, slope=lm.psr.mdd.pe,color="cornflowerblue")
      + geom_abline(intercept=lm.pom.sad.intercept, slope=lm.pom.sad.pe,color="coral")
      + geom_abline(intercept=lm.pom.tin.intercept, slope=lm.pom.tin.pe,color="chartreuse3")
      + geom_abline(intercept=lm.pom.mdd.intercept, slope=lm.pom.mdd.pe,color="cornflowerblue")
  
#### Figure 3 -- ROC curves
    
    library(pROC)
  
  
    # Check integrity of response variables
    describeBy(good.treatment.wcst$PSR,good.treatment.wcst$respondent25)
    describeBy(good.treatment.wcst$PSR,good.treatment.wcst$respondent50)
    describeBy(good.treatment.wcst$PSR,good.treatment.wcst$respondent75)

  
  ## Full sample
  
    # 75%
  
    roc.all.r75.adjPE = roc(wcst=good.treatment.wcst, good.treatment.wcst$respondent75, good.treatment.wcst$adjPE,auc = T,plot=F,adjPErcent=T)
    roc.all.r75.adjPE
    set.seed(2015)
    ci.auc(roc.all.r75.adjPE, method="b")
    coords(roc.all.r75.adjPE, "best")
    ci.thresholds(roc.all.r75.adjPE,thresholds = "best")
    has.partial.auc(roc.all.r75.adjPE)
  
    # 50%
    
    roc.all.r50.adjPE = roc(wcst=good.treatment.wcst, good.treatment.wcst$respondent50, good.treatment.wcst$adjPE,auc = T,plot=F,adjPErcent=T)
    roc.all.r50.adjPE
    set.seed(2015)
    ci.auc(roc.all.r50.adjPE, method="b")
    coords(roc.all.r50.adjPE, "best")
    ci.thresholds(roc.all.r50.adjPE,thresholds = "best")
    has.partial.auc(roc.all.r50.adjPE)
  
    # 25%
    
    roc.all.r25.adjPE = roc(wcst=good.treatment.wcst, good.treatment.wcst$respondent25, good.treatment.wcst$adjPE,auc = T,plot=F,adjPErcent=T)
    roc.all.r25.adjPE
    set.seed(2015)
    ci.auc(roc.all.r25.adjPE, method="b")
    coords(roc.all.r25.adjPE, "best")
    ci.thresholds(roc.all.r25.adjPE,thresholds = "best")
    has.partial.auc(roc.all.r25.adjPE)
  
    # Full sample plot
    
    plot(roc.all.r75.adjPE,col="gold",lwd=4,grid=T,add=F)
    plot(roc.all.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.all.r25.adjPE,col="chocolate3",add=T,lwd=4)

  ## SAD sample
    
    # 75%
  
    roc.sad.r75.adjPE = roc(wcst=good.sad.treatment.wcst, good.sad.treatment.wcst$respondent75, good.sad.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.sad.r75.adjPE
    set.seed(2015)
    ci.auc(roc.sad.r75.adjPE, method="b")
    coords(roc.sad.r75.adjPE, "best")
    ci.thresholds(roc.sad.r75.adjPE,thresholds = "best")
    has.partial.auc(roc.sad.r75.adjPE)
  
    # 50%
    
    roc.sad.r50.adjPE = roc(wcst=good.sad.treatment.wcst, good.sad.treatment.wcst$respondent50, good.sad.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.sad.r50.adjPE
    set.seed(2015)
    ci.auc(roc.sad.r50.adjPE, method="b")
    coords(roc.sad.r50.adjPE, "best")
    ci.thresholds(roc.sad.r50.adjPE,thresholds = "best")
    has.partial.auc(roc.sad.r50.adjPE)
  
    # 25%
    
    roc.sad.r25.adjPE = roc(wcst=good.sad.treatment.wcst, good.sad.treatment.wcst$respondent25, good.sad.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.sad.r25.adjPE
    set.seed(2015)
    ci.auc(roc.sad.r25.adjPE, method="b")
    coords(roc.sad.r25.adjPE, "best")
    ci.thresholds(roc.sad.r25.adjPE,thresholds = "best")
    has.partial.auc(roc.sad.r25.adjPE)
  
    # SAD plot
    
    plot(roc.sad.r75.adjPE,col="gold",lwd=4,grid=T,add=F)
    plot(roc.sad.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.sad.r25.adjPE,col="chocolate3",add=T,lwd=4)
    
  ## Tinnitus sample
    
    # 75%
    
    roc.tin.r75.adjPE = roc(wcst=good.tin.treatment.wcst, good.tin.treatment.wcst$respondent75, good.tin.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.tin.r75.adjPE
    set.seed(2015)
    ci.auc(roc.tin.r75.adjPE, method="b")
    coords(roc.tin.r75.adjPE, "best")
    ci.thresholds(roc.tin.r75.adjPE,thresholds = "best")
    has.partial.auc(roc.tin.r75.adjPE)
    
    # 50%
    
    roc.tin.r50.adjPE = roc(wcst=good.tin.treatment.wcst, good.tin.treatment.wcst$respondent50, good.tin.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.tin.r50.adjPE
    set.seed(2015)
    ci.auc(roc.tin.r50.adjPE, method="b")
    coords(roc.tin.r50.adjPE, "best")
    ci.thresholds(roc.tin.r50.adjPE,thresholds = "best")
    has.partial.auc(roc.tin.r50.adjPE)
    
    # 25%
    
    roc.tin.r25.adjPE = roc(wcst=good.tin.treatment.wcst, good.tin.treatment.wcst$respondent25, good.tin.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.tin.r25.adjPE
    set.seed(2015)
    ci.auc(roc.tin.r25.adjPE, method="b")
    coords(roc.tin.r25.adjPE, "best")
    ci.thresholds(roc.tin.r25.adjPE,thresholds = "best")
    has.partial.auc(roc.tin.r25.adjPE)
    
    # TIN plot
    
    plot(roc.tin.r75.adjPE,col="gold",lwd=4,grid=T,add=F)
    plot(roc.tin.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.tin.r25.adjPE,col="chocolate3",add=T,lwd=4)
    
  ## MDD sample
    
    # 75%
    
    roc.mdd.r75.adjPE = roc(wcst=good.mdd.treatment.wcst, good.mdd.treatment.wcst$respondent75, good.mdd.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.mdd.r75.adjPE
    set.seed(2015)
    ci.auc(roc.mdd.r75.adjPE, method="b")
    coords(roc.mdd.r75.adjPE, "best")
    ci.thresholds(roc.mdd.r75.adjPE,thresholds = "best")
    has.partial.auc(roc.mdd.r75.adjPE)
    
    # 50%
    
    roc.mdd.r50.adjPE = roc(wcst=good.mdd.treatment.wcst, good.mdd.treatment.wcst$respondent50, good.mdd.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.mdd.r50.adjPE
    set.seed(2015)
    ci.auc(roc.mdd.r50.adjPE, method="b")
    coords(roc.mdd.r50.adjPE, "best")
    ci.thresholds(roc.mdd.r50.adjPE,thresholds = "best")
    has.partial.auc(roc.mdd.r50.adjPE)
    
    # 25%
    
    roc.mdd.r25.adjPE = roc(wcst=good.mdd.treatment.wcst, good.mdd.treatment.wcst$respondent25, good.mdd.treatment.wcst$adjPE, auc = T,plot=F,adjPErcent=T)
    roc.mdd.r25.adjPE
    set.seed(2015)
    ci.auc(roc.mdd.r25.adjPE, method="b")
    coords(roc.mdd.r25.adjPE, "best")
    ci.thresholds(roc.mdd.r25.adjPE,thresholds = "best")
    has.partial.auc(roc.mdd.r25.adjPE)
    
    # MDD plot
    
    plot(roc.mdd.r75.adjPE,col="gold",lwd=4,grid=T,add=F)
    plot(roc.mdd.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.mdd.r25.adjPE,col="chocolate3",add=T,lwd=4)

    
  ## Plot all (export as 14*14 PDF)

    par(mfrow=c(2,2),bg="white")
    plot(roc.all.r75.adjPE,col="gold",lwd=4,grid=T,add=F,main="Full sample")
    plot(roc.all.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.all.r25.adjPE,col="chocolate3",add=T,lwd=4)

    plot(roc.sad.r75.adjPE,col="gold",lwd=4,grid=T,add=F,main="Social anxiety disorder sample")
    plot(roc.sad.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.sad.r25.adjPE,col="chocolate3",add=T,lwd=4)
    
    plot(roc.tin.r75.adjPE,col="gold",lwd=4,grid=T,add=F,main="Tinnitus sample")
    plot(roc.tin.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.tin.r25.adjPE,col="chocolate3",add=T,lwd=4)
    
    plot(roc.mdd.r75.adjPE,col="gold",lwd=4,grid=T,add=F,main="Depression sample")
    plot(roc.mdd.r50.adjPE,col="cornsilk3",add=T,lwd=4)
    plot(roc.mdd.r25.adjPE,col="chocolate3",add=T,lwd=4)
    
    par(old.par)