# R code for lecture V_CAN4_08 ------------------------------------------------ 

# packages --------------------------------------------------------------------

# you may need to install them first using (in RStudio) Tools>Install Packages
library(psych)      # for diverse helper functions
library(robustbase) # for robust regression


# figure speed-accuracy trade off ----
x=seq(-3,3,by=.01)
par(mar=c(5,5,4,2))
plot(x,pnorm(x),type="l",xlab="Response Time",ylab="Accuracy",cex.lab=1.5,lwd=2,axes=F)
axis(1)
axis(2)

# figure RT variability ----

RT=NULL
for (i in 1:105) {
  rt=scale(rlnorm(200,0,.5))
  rt=rt*sample(rnorm(500,25,5),1)+sample(rnorm(500,500,100),1)
  RT=cbind(RT,rt)
}
mRT =apply(RT,2,mean)
sdRT=apply(RT,2,sd)
cvRT=sdRT/mRT
summary(lm(sdRT~mRT))

plot(RT[,1],type="l",xlab="Trial",ylab="Response Time",xlim=c(1,225),cex.lab=1.5,col=grey(.75))
lines(par("usr")[1:2],rep(mean(RT[,1]),2),col=2,lwd=2)
lines(par("usr")[1:2],rep(mean(RT[,1])+sd(RT[,1]),2),col=2,lwd=2,lty=3)
lines(par("usr")[1:2],rep(mean(RT[,1])-sd(RT[,1]),2),col=2,lwd=2,lty=3)
rect(205,par("usr")[3],par("usr")[2],par("usr")[4],col=0,border=NA)
text(rep(220,3),c(mean(RT[,1])+sd(RT[,1]),mean(RT[,1]),mean(RT[,1])-sd(RT[,1])),c("+1SD","M","-1SD"))
box()
par(mar=c(5,4,4,2))

# data ----
summary_stroop_2014_12_09 <- read.delim2("~/Documents/R/CAN4/summary_stroop_2014_12_09.txt")
stroop=summary_stroop_2014_12_09
# the data contain values for different RT-derives measures:
# standard (in sec)
# - simple RT 
# - coeficient of variation
# ex-gauss fit parameters (in msec, need to be rescaled later on)
# - mu
# - sigma
# - tau
# diffusion modeling parameters (in sec/arbitrary units a.u.)
# - v (drift rate)
# - a (boundary separation)
# - Ter (indecision time)

# basic stats ----
pc_c= stroop$Pc_C              # % correct congruent stroop trials
pc_i= stroop$Pc_I              # same for incongruent trials 
m_c = stroop$MRT_C             # mean RT
m_i = stroop$MRT_I             
cv_c= sqrt(stroop$VRT_C)/m_c   # coefficient of variation = sd RT / mean RT (VRT_x is variance, hence sqrt)
cv_i= sqrt(stroop$VRT_I)/m_i
performance=data.frame(pc_c,pc_i,m_c,m_i,cv_c,cv_i)

# data overview
# we want to have a 2 row and two column plot
par(mfrow=c(2,2))
# we want - of course - box plots
b0=boxplot(performance,ylim=c(0,2),main="Boxplot")
# returns Mahalanobis distances (MD) to detect multivariate outliers 
# and plots them against their associated Chi-Square quantiles
o0=outlier(performance,ylim=c(0,100))
# is Mahalanobis distance greater that its associated 99.9 quantile 
# of the Chi-Square distribution with as much df as there are vars?
out=which(o0>qchisq(.999,ncol(performance))) 
# if so, leave out respective cases via v[-out] (retains only those 
# cases that are not in 'out' vector)
pc_c= stroop$Pc_C[-out] 
pc_i= stroop$Pc_I[-out]
m_c = stroop$MRT_C[-out]
m_i = stroop$MRT_I[-out]
cv_c= sqrt(stroop$VRT_C[-out])/m_c
cv_i= sqrt(stroop$VRT_I[-out])/m_i
performance=data.frame(pc_c,pc_i,m_c,m_i,cv_c,cv_i)
# plot what the data look line without outliers 
b1=boxplot(performance,ylim=c(0,2),main="Boxplot")
o1=outlier(performance,ylim=c(0,100))
out1=which(o1>qchisq(.999,ncol(performance)))
# looks way better, still not optimal, but this ist for exercise and 
# visualization purposes only 
par(mfrow=c(1,1)) # reset number of plots

# some more basic stats ----

# exgauss fit
mu_c=stroop$R_mu_C[-out]
mu_i=stroop$R_mu_I[-out]
sig_c=stroop$R_sigma_C[-out]
sig_i=stroop$R_sigma_I[-out]
tau_c=stroop$R_tau_C[-out]
tau_i=stroop$R_tau_I[-out]
# EZ drift diffusion parameters
v_c=stroop$v_C[-out]
v_i=stroop$v_I[-out]
a_c=stroop$a_C[-out]
a_i=stroop$a_I[-out]
Ter_c=stroop$Ter_C[-out]
Ter_i=stroop$Ter_I[-out]

performance2=data.frame(pc_c,pc_i,m_c,m_i,cv_c,cv_i,mu_c,mu_i,sig_c,sig_i,tau_c,tau_i,v_c,v_i,a_c,a_i,Ter_c,Ter_i)
performance2[,7:12]=performance2[,7:12]/1000

# relationship RT/Percent Correct ----
par(mar=c(5,5,4,2))
plot(m_c~pc_c,performance2,col=grey(.75),xlim=c(.8,1),ylim=c(.4,1),las=1,xlab="Percent Correct",ylab="Response Time",cex.lab=1.5)
abline(lmrob(m_c~pc_c,performance2),col=grey(.75),lwd=2) 
points(m_i~pc_i,performance2,col=grey(.50),pch=19)
abline(lmrob(m_i~pc_i,performance2),col=grey(.50),lwd=2)
legend.text=c("Robust Regression",
              paste("Con (beta = ",round(lmrob(m_c~pc_c,data.frame(scale(performance)))$coefficients[2],2),")",sep=""),
              paste("Inc (beta = ",round(lmrob(m_i~pc_i,data.frame(scale(performance)))$coefficients[2],2),")",sep=""))
legend("topleft",inset=.01,legend=legend.text,lty=c(NA,1,1),lwd=c(NA,2,2),col=c(NA,grey(.75),grey(.50)),bty="n",cex=1.15)
par(mar=c(5,4,4,2))

# plot correlations with Percent Correct --------------------

# preliminaries
correlations=NULL
for (i in seq(3,17,2)) {
  correlations=rbind(correlations,c(cor(performance[,1],performance[,i],method="spearman"),cor(performance[,2],performance[,i+1],method="spearman")))
}
colnames(correlations)=c("Congruent","Incongruent")
rownames(correlations)=c("Mean RT","CV RT","mu","sigma","tau","v","a","Ter")

# CI correlations
# function for calculating 95% CI of a given correlation at a given sample size 
# what we need to do is to transform the correlation to z scale via the fisher z
# transformation using psych::fisherz, calculate 95% CI via qnorm((.975) times the
# approx. standard error sqrt(1/(n-3))and then backtransform the result via 
# psych:: fisherz2r
cor.ci <- function(cor,n,ci.only=F) {
  fz  = fisherz(cor)
  dev = qnorm(.975)*sqrt(1/(n-3))
  upr = fisherz2r(fz+dev)
  lwr = fisherz2r(fz-dev)
  if (ci.only==F) {
    out = c(lwr,cor,upr)
  } else {
    out = c(lwr,upr)
  }
  return(out)
}

par(mfrow=c(3,3))   # plot 3x3 panels
par(mar=c(4,4,1,2)) # with smaller than the default margins
# mean RT (comment: las=1 plots alls axis tick labels horizontal)
plot(m_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(.4,1),las=1,xlab="Percent Correct",ylab="Response Time")
abline(lmrob(m_c~pc_c,performance2),col=grey(.75)) 
# robst regression via robustbase::lmrob to avoid influence of possible unidentified outliers
# abline uses the coefficients (i.e. intercept and slope) returned by lmrob to plot the regression line 
points(m_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(m_i~pc_i,performance2),col=grey(.50))
# CV RT
plot(cv_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(.15,1),las=1,xlab="Percent Correct",ylab="Coefficient of Variation")
abline(lmrob(cv_c~pc_c,performance2),col=grey(.75))
points(cv_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(cv_i~pc_i,performance2),col=grey(.50))
# legend for whole plot
legend("topleft",inset=.01,bty="n",legend=c("Congruent","Incongruent"),lty=1,col=c(grey(.75),grey(.50)),title="Robust Fit")
# correlation bar plot
# we could use barplot instead, but this gives us limited control
plot(c(min(correlations),max(correlations)),c(1,8),ylim=c(.5,8.5),type="n",axes=F,las=1,ylab="",xlab="Spearman Correlation with Percent Correct",sub="Shading indicates 95% CI of zero correlation")
ci.0=cor.ci(0,nrow(performance),ci.only=T)
rect(ci.0[1],par("usr")[3],ci.0[2],par("usr")[4],col="#FFCCCC",border=NA)
axis(1)
lines(rep(0,2),par("usr")[3:4],lty=1)
for (i in 1:8) {
  rect(0,9-i+.25,correlations[i,1],9-i,col=grey(.75))
  rect(0,9-i,correlations[i,2],9-i-.25,col=grey(.50))
}
axis(2,at=8:1,c("RT","CV RT","mu","sigma","tau","v","a","Ter"),las=1,tick=F)

# mu
plot(mu_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(.3,.8),las=1,xlab="Percent Correct",ylab="Ex-Gauss: Mu")
abline(lmrob(mu_c~pc_c,performance2),col=grey(.75))
points(mu_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(mu_i~pc_i,performance2),col=grey(.50))
# sigma
plot(sig_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(.0,.2),las=1,xlab="Percent Correct",ylab="Ex-Gauss: Sigma")
abline(lmrob(sig_c~pc_c,performance2),col=grey(.75))
points(sig_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(sig_i~pc_i,performance2),col=grey(.50))
# tau
plot(tau_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(.0,.8),las=1,xlab="Percent Correct",ylab="Ex-Gauss: Tau")
abline(lmrob(tau_c~pc_c,performance2),col=grey(.75))
points(tau_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(tau_i~pc_i,performance2),col=grey(.50))
# v
plot(v_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(.0,.5),las=1,xlab="Percent Correct",ylab="Diffusion: v")
abline(lmrob(v_c~pc_c,performance2),col=grey(.75))
points(v_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(v_i~pc_i,performance2),col=grey(.50))
# a
plot(a_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(.0,.5),las=1,xlab="Percent Correct",ylab="Diffusion: a")
abline(lmrob(a_c~pc_c,performance2),col=grey(.75))
points(a_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(a_i~pc_i,performance2),col=grey(.50))
# Ter
plot(Ter_c~pc_c,performance2,col=grey(.75),cex=.5,xlim=c(.8,1),ylim=c(-.1,.5),las=1,xlab="Percent Correct",ylab="Diffusion: Ter")
abline(lmrob(Ter_c~pc_c,performance2),col=grey(.75))
points(Ter_i~pc_i,performance2,col=grey(.50),cex=.5,pch=19)
abline(lmrob(Ter_i~pc_i,performance2),col=grey(.50))
par(mfrow=c(1,1))

# differences between congruent and incongruent using div RT derived measures ----

# calculates means (ms) and SEM (es) for all variables of interest
ms2=colMeans(performance2)
ss2=apply(performance2,2,sd)
es2=ss2/sqrt(nrow(performance2))

# function for bootstrapped difference
# why bootstrap? simply because (1) assumptions of parametric procedures 
# may not be met and (2) Mann Whitney U test for nonparametric test sounds boring 
boot.difference <- function(x1,x2,replicates=1000) {
  set.seed(242)
  # data
  df=data.frame(x1,x2)
  colnames(df)=c("x1","x2")
  
  # Bootstrap 95% CI for correlation coefficients 
  library(boot)
  
  # function to obtain correlation coefficient 
  ds = function(data, indices) {
    d  = data[indices,] # allows boot to select sample 
    diff  = mean(d$x2)-mean(d$x1)
    return(diff) 
  } 
  
  # bootstrapping with 1000 replications
  results  <- boot(data=df, statistic=ds, R=replicates, parallel="multicore", ncpus=4)
  
  # get 95% confidence intervals 
  bci=boot.ci(results, type="bca", index=1)
  CI.lo=bci$bca[4]
  CI.hi=bci$bca[5]
  
  # return result
  bootMean=mean(results$t)
  bootSE=sd(results$t)
  z=bootMean/bootSE
  p=1-pnorm(abs(z))
  s=data.frame(bootMean,bootSE,z,p,CI.lo,CI.hi)
  return(s)
  
}
b.pc  = boot.difference(performance$pc_c,performance$pc_i)
b.m   = boot.difference(performance$m_c,performance$m_i)
b.cv  = boot.difference(performance$cv_c,performance$cv_i)
b.mu  = boot.difference(performance2$mu_c,performance2$mu_i)
b.sig = boot.difference(performance2$sig_c,performance2$sig_i)
b.tau = boot.difference(performance2$tau_c,performance2$tau_i)
b.v   = boot.difference(performance2$v_c,performance2$v_i)
b.a   = boot.difference(performance2$a_c,performance2$a_i)
b.Ter = boot.difference(performance2$Ter_c,performance2$Ter_i)
boot.results=rbind(b.pc,b.m,b.cv,b.mu,b.sig,b.tau,b.v,b.a,b.Ter)
rownames(boot.results)=c("Correct","Mean RT","CV RT","mu","sigma","tau","v","a","Ter")


par(mar=c(5,5,3,5))
plot(c(1,9),c(0,1),type="n",axes=F,ylab="Correct Responses [Frequency]",xlab="",xlim=c(.5,9.5),cex.lab=1.5)
axis(1,at=1:9,c("Correct","RT","CV RT","mu","sigma","tau","v","a","Ter"),tick=F)
axis(2,las=1)
axis(4,las=1) # plots another axis to the right of the plot
mtext(side = 4, line = 3, "Response Time Measures [s/a.u.]",cex=1.5) # mtext lets you plot in the margin
for (i in seq(1,17,2)) {
  j=ceiling(i/2)
  rect(j-.25,0,j,ms2[i],col=grey(.75))
  lines(rep(j-.125,2),c(ms2[i]-es2[i],ms2[i]+es2[i]))
  rect(j,0,j+.25,ms2[i+1],col=grey(.50))
  lines(rep(j+.125,2),c(ms2[i+1]-es2[i+1],ms2[i+1]+es2[i+1]))
  p=""
  if (boot.results$p[j]<.05) {
    p="*"
  }
  if (boot.results$p[j]<.01) {
    p="**"
  }
  if (boot.results$p[j]<.001) {
    p="***"
  }
  text(j,1,p,cex=1.5)
}
lines(c(2,3),rep(-.2,2),lwd=2,xpd=T)      # xpd=T lets you plot the line in the margin
mtext("Standard",side=1,line=3.25,at=2.5) 
lines(c(4,6),rep(-.2,2),lwd=2,xpd=T)
mtext("Ex-Gauss Fitting",side=1,line=3.25,at=5)
lines(c(7,9),rep(-.2,2),lwd=2,xpd=T)
mtext("Diffusion Modeling",side=1,line=3.25,at=8)

# "Exkurs im Exkurs": Illustration of behavior of linear vs. robust regression in the presence of possible outliers 
# tau_i 
df=data.frame(performance2$tau_i,performance2$pc_i)
colnames(df)=c("tau_i","pc_i")
out.df=which(outlier(df,plot=F)>10) # retrieves possible outliers
plot(tau_i~pc_i,performance2,col=grey(.75),las=1,xlab="Percent Correct",ylab="Ex-Gauss: Tau",cex.lab=1.5)
points(tau_i~pc_i,df[out.df,],col=1,pch=19) # marks possible outliers
abline(lm(tau_i~pc_i,df),col=grey(.75),lwd=2) # linear regression with possible outliers
abline(lmrob(tau_i~pc_i,df),col=2,lwd=2) # robust regression with possible outliers
abline(lm(tau_i~pc_i,df[-out.df,]),col=grey(.75),lwd=2,lty=2) # linear regression without possible outliers
abline(lmrob(tau_i~pc_i,df[-out.df,]),col=2,lwd=2,lty=2) # robust regression without possible outliers