Reproduce results of rmcorr paper

Jon Bakdash and Laura Marusich

February 6, 2017

###1. Figure 1: rmcorr and reg plot

###2. Figure 2: rmcorr vs OLS reg

###3. Figure 3: rmcorr w/data transformations

###4. Figure 4: Power curves 

## png 
##   2
## null device 
##           1

###5. Brain volume and age rmcorr and simple reg/cor results and Figure 5 #### a) rmcorr and simple reg results

## Warning in rmcorr(participant = Participant, measure1 = Age, measure2 =
## Volume, : 'Participant' coerced into a factor
## 
## Repeated measures correlation
## 
## r
## -0.7044077
## 
## degrees of freedom
## 71
## 
## p-value
## 3.561007e-12
## 
## 95% confidence interval
## -0.8053581 -0.5637514
## 
##  Pearson's product-moment correlation
## 
## data:  avgRaz2005$Age and avgRaz2005$Volume
## t = -3.4912, df = 70, p-value = 0.000837
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5662456 -0.1684542
## sample estimates:
##        cor 
## -0.3850943
## 
## Call:
## lm(formula = Volume ~ Age, data = raz2005)
## 
## Coefficients:
## (Intercept)          Age  
##    151.9068      -0.3399
## 
##  Pearson's product-moment correlation
## 
## data:  raz2005$Age and raz2005$Volume
## t = -5.165, df = 142, p-value = 7.984e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5269809 -0.2503991
## sample estimates:
##        cor 
## -0.3976861

b) Figure 5: Brain vol + age, rmcorr and simple reg/cor

## png 
##   2
## null device 
##           1

###6. Visual search rmcorr and simple reg/cor results and Figure 6 #### a) rmcorr and simple reg results

## Warning in rmcorr(participant = sub, measure1 = rt, measure2 = acc, dataset =
## gilden2010): 'sub' coerced into a factor
## 
## Repeated measures correlation
## 
## r
## -0.406097
## 
## degrees of freedom
## 32
## 
## p-value
## 0.01716871
## 
## 95% confidence interval
## -0.6611673 -0.06687244
## 
## Call:
## lm(formula = acc ~ rt, data = gildenMeans)
## 
## Coefficients:
## (Intercept)           rt  
##      0.8132       0.1777
## 
##  Pearson's product-moment correlation
## 
## data:  gildenMeans$rt and gildenMeans$acc
## t = 2.1966, df = 9, p-value = 0.05565
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.01409542  0.87910346
## sample estimates:
##       cor 
## 0.5907749
## 
## Call:
## lm(formula = acc ~ rt, data = gilden2010)
## 
## Coefficients:
## (Intercept)           rt  
##      0.8612       0.1111
## 
##  Pearson's product-moment correlation
## 
## data:  gilden2010$rt and gilden2010$acc
## t = 2.6401, df = 42, p-value = 0.01158
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.09053513 0.60625185
## sample estimates:
##       cor 
## 0.3772751

b) Figure 6: Visual search rmcorr and simple reg/cor

## png 
##   2
## null device 
##           1

##Appendix A: Rmcorr minus standard corr power

## png 
##   2
## null device 
##           1
  1. Rmcorr and multi-level model with Raz et al. 2005 data
brainvolage.rmc <- rmcorr(participant = Participant, measure1 = Age, measure2 = Volume, dataset = raz2005)
## Warning in rmcorr(participant = Participant, measure1 = Age, measure2 =
## Volume, : 'Participant' coerced into a factor
#Null multi-level model: Random intercept and fixed slope
null.vol <- lmer(Volume ~ Age + (1 | Participant), data = raz2005, REML = FALSE)

#Model fit
null.vol
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: Volume ~ Age + (1 | Participant)
##    Data: raz2005
##       AIC       BIC    logLik  deviance  df.resid 
##  977.4705  989.3497 -484.7352  969.4705       140 
## Random effects:
##  Groups      Name        Std.Dev.
##  Participant (Intercept) 11.287  
##  Residual                 3.024  
## Number of obs: 144, groups:  Participant, 72
## Fixed Effects:
## (Intercept)          Age  
##    163.7090      -0.5533
#Parameter Confidence Intervals
confint(null.vol)
## Computing profile confidence intervals ...
##                   2.5 %      97.5 %
## .sig01        9.5208850  13.6036847
## .sigma        2.5713387   3.6236647
## (Intercept) 155.1479333 172.3796868
## Age          -0.7016833  -0.4056009
#Model fitted values and confidence intervals for each participant (L1 effects)
set.seed(9999)
L1.predict.raz <- predictInterval(null.vol, newdata = raz2005, n.sims = 1000)
L1.predict.raz <- cbind(raz2005$Participant, L1.predict.raz)

apatheme =  theme_bw() +
            theme(panel.grid.major=element_blank(),
                  panel.grid.minor=element_blank(),
                  panel.border=element_blank(),
                  axis.line=element_line(),
                  legend.position="none",
                  #text=element_text(family='Times'), #Need- https://github.com/wch/extrafont + Ghostscript
                  #legend.title=element_blank(),
                  axis.line.x = element_line(color="black", size = 0.9),
                  axis.line.y = element_line(color="black", size = 0.9),
                  axis.text.x = element_text(size = 12),
                  axis.text.y = element_text(size = 12), 
                  axis.title.x = element_text(size = 12),
                  axis.title.y = element_text(size = 12))

#Create custom color palette 
Blues<-brewer.pal(9,"Blues")

ggplot(raz2005, aes(x = Age, y = Volume, group = Participant, color = Participant)) +
  geom_line(aes(y = predict(null.vol)), linetype = 2) + 
  geom_line(aes(y = brainvolage.rmc$model$fitted.values), linetype = 1) + 
  geom_ribbon(aes(ymin = L1.predict.raz$lwr,
                  ymax = L1.predict.raz$upr,
                  group = L1.predict.raz$`raz2005$Participant`,
                  linetype = NA), alpha = 0.07) +
  apatheme + 
  labs(title = "Rmcorr and Random Intercept Multi-Level Model:\n Raz et al. 2005 Data", x = "Age", 
       y = expression(Cerebellar~Hemisphere~Volume~(cm^{3}))) +
  geom_point(aes(colour = Participant)) +
  scale_colour_gradientn(colours=Blues) + theme(plot.title = element_text(hjust = 0.5)) + 
  geom_abline(intercept = fixef(null.vol)[1], slope = fixef(null.vol)[2], colour = "black", size = 1, linetype = 2)

ggsave(file = "plots/AppendixC_Figure1.eps", width = 5.70 , height = 5.73, dpi = 300)
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page
  1. Rmcorr and multi-level model with Gilden et al. 2010 data
vissearch.rmc <- rmcorr(participant = sub, measure1 = rt, measure2 = acc, dataset = gilden2010)
## Warning in rmcorr(participant = sub, measure1 = rt, measure2 = acc, dataset =
## gilden2010): 'sub' coerced into a factor
null.vis <- lmer(acc ~ rt + (1 | sub), data = gilden2010, REML = FALSE)

#Model 1: Random intercept + random slope for RT
model1.vis <- lmer(acc ~ rt + (1 + rt | sub), data = gilden2010, REML = FALSE)
## boundary (singular) fit: see ?isSingular
#Model Comparison
#a) Chi-Square
anova(null.vis, model1.vis)
## Data: gilden2010
## Models:
## null.vis: acc ~ rt + (1 | sub)
## model1.vis: acc ~ rt + (1 + rt | sub)
##            npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## null.vis      4 -197.41 -190.27 102.70  -205.41                     
## model1.vis    6 -193.46 -182.75 102.73  -205.46 0.0497  2     0.9754
# #b) Evidence ratio using AIC
# Models.vis<-list()
# Models.vis<-c(null.vis, model1.vis)
# 
# ModelTable2<-AICcmodavg::aictab(Models.vis, modnames = c("null", "Model 1"))
# ModelTable2
# evidence(ModelTable2)

  #Estimating CIs: Convergence problems with model 1
  confint(model1.vis)
## Computing profile confidence intervals ...
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep
## Warning in FUN(X[[i]], ...): non-monotonic profile for .sig01
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep

## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): unexpected decrease in
## profile: using minstep
## Warning in FUN(X[[i]], ...): non-monotonic profile for .sig03
## Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit
## for .sig01: falling back to linear interpolation
## Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit
## for .sig03: falling back to linear interpolation
##                   2.5 %     97.5 %
## .sig01       0.02197527 0.11293808
## .sig02      -1.00000000 1.00000000
## .sig03       0.00000000        Inf
## .sigma       0.01303108 0.02157782
## (Intercept)  0.91458725 1.06378202
## rt          -0.15425292 0.03364490
  warnings()
  
  set.seed(9999)
  predictInterval(model1.vis, newdata = gilden2010, n.sims = 1000)
##          fit       upr       lwr
## 1  0.8663920 0.8960908 0.8352832
## 2  0.8683327 0.8992494 0.8385416
## 3  0.8707389 0.8997923 0.8393785
## 4  0.8725921 0.9043416 0.8411908
## 5  0.9429413 0.9732029 0.9159473
## 6  0.9432319 0.9722811 0.9162444
## 7  0.9500212 0.9787765 0.9220920
## 8  0.9521220 0.9792167 0.9237837
## 9  0.9481541 0.9752953 0.9207015
## 10 0.9557798 0.9844996 0.9267692
## 11 0.9559257 0.9853144 0.9268080
## 12 0.9596370 0.9887583 0.9285575
## 13 0.9443758 0.9730895 0.9144351
## 14 0.9430301 0.9705867 0.9121712
## 15 0.9481910 0.9777425 0.9189827
## 16 0.9500669 0.9785099 0.9180516
## 17 0.9546490 0.9828561 0.9234057
## 18 0.9511326 0.9826628 0.9211205
## 19 0.9583308 0.9867742 0.9280774
## 20 0.9601212 0.9901520 0.9312060
## 21 0.9227967 0.9488316 0.8935664
## 22 0.9244902 0.9553428 0.8959148
## 23 0.9277636 0.9569939 0.8986043
## 24 0.9300477 0.9599092 0.9000850
## 25 0.9625122 0.9928161 0.9347764
## 26 0.9680937 0.9978023 0.9390449
## 27 0.9706817 1.0010774 0.9416283
## 28 0.9749069 1.0033890 0.9446229
## 29 0.9337347 0.9613853 0.9051130
## 30 0.9427611 0.9688831 0.9133651
## 31 0.9495472 0.9767039 0.9189270
## 32 0.9506709 0.9813415 0.9226078
## 33 0.9092402 0.9374052 0.8801574
## 34 0.9081569 0.9380428 0.8783812
## 35 0.9101271 0.9393227 0.8784178
## 36 0.9130451 0.9411345 0.8840515
## 37 0.9530667 0.9841683 0.9233226
## 38 0.9633384 0.9921144 0.9363187
## 39 0.9616431 0.9909585 0.9335434
## 40 0.9655959 0.9935503 0.9362402
## 41 0.9739775 1.0029519 0.9440092
## 42 0.9725509 1.0019449 0.9447089
## 43 0.9767107 1.0060772 0.9455765
## 44 0.9756523 1.0049782 0.9452859
  warnings()

#Estimating and graphing null model 
#Parameter Confidence Intervals
confint(null.vis)
## Computing profile confidence intervals ...
##                   2.5 %     97.5 %
## .sig01       0.02138120 0.05796124
## .sigma       0.01294851 0.02139924
## (Intercept)  0.91815424 1.05718722
## rt          -0.15404840 0.02955915
#Model fitted values and confidence intervals for each participant (L1 effects)
set.seed(9999)
L1.predict.gilden <- predictInterval(null.vis, newdata = gilden2010, n.sims = 1000)
L1.predict.gilden <- cbind(gilden2010$sub, L1.predict.gilden)

apatheme =  theme_bw() +
            theme(panel.grid.major=element_blank(),
                  panel.grid.minor=element_blank(),
                  panel.border=element_blank(),
                  axis.line=element_line(),
                  legend.position="none",
                  #text=element_text(family='Times'), #Need- https://github.com/wch/extrafont + Ghostscript
                  #legend.title=element_blank(),
                  axis.line.x = element_line(color="black", size = 0.9),
                  axis.line.y = element_line(color="black", size = 0.9),
                  axis.text.x = element_text(size = 12),
                  axis.text.y = element_text(size = 12), 
                  axis.title.x = element_text(size = 12),
                  axis.title.y = element_text(size = 12))

#Create custom color palette 
Colors12<-brewer.pal(12,"Paired")

ggplot(gilden2010, aes(x = rt, y = acc, group = sub, color = sub)) +
  geom_line(aes(y = predict(null.vis)), linetype = 2) + 
  geom_line(aes(y = vissearch.rmc$model$fitted.values), linetype = 1) + 
  geom_ribbon(aes(ymin = L1.predict.gilden$lwr,
                  ymax = L1.predict.gilden$upr,
                  group = L1.predict.gilden$`gilden2010$sub`,
                  linetype = NA),
                  alpha = 0.07) +
  apatheme + theme(plot.title = element_text(hjust = 0.5)) + 
  labs(title = "Rmcorr and Random Intercept Multi-Level Model:\n Gilden et al. 2010 Data", x = "Response Time (Seconds)", y = "Accuracy") +
  geom_point(aes(colour = sub)) +
  scale_colour_gradientn(colours=Colors12) + 
  geom_abline(intercept = fixef(null.vis)[1], slope = fixef(null.vis)[2], colour = "black", size = 1, linetype = 2) +
  
  scale_y_continuous(breaks=seq(0.80, 1.0, 0.05)) +
  scale_x_continuous(breaks=seq(0.50, 0.9, 0.1)) 

ggsave(file = "plots/AppendixC_Figure2.eps", width = 5.70 , height = 5.73, dpi = 300)
## Warning in grid.Call.graphics(C_polygon, x$x, x$y, index): semi-transparency is
## not supported on this device: reported only once per page