library(lsr)            # compute effect sizes in ttests
library(R.matlab)
library(plyr)           # to summarize data and required before dplyr
library(dplyr)          # for renaming columns 
library(reshape)        # melt and cast functions
library(BayesFactor)
library(knitr)

source(paste(fullPath, "/functions/as.data.frame.list.R", sep=""))
# Load all "preprocessed" data
load(file = paste(fullPath,'/data/modelresults.Rda', sep=""))
load(file = paste(fullPath,'/data/modelresults_long.Rda', sep=""))
load(file = paste(fullPath,'/data/congruencyIndex.Rda', sep=""))
load(file = paste(fullPath,'/data/resultsRH_perSubj.Rda', sep=""))
load(file = paste(fullPath,'/data/resultsRH_perTrial.Rda', sep=""))

1 Analysis of Sequence Effects

1.1 r parameter (number of trials in line with RH)

# Do an aov first
modelResults.summary <- aggregate(r_param ~ conditionOrder * taskOrder * condition * ID, modelResults.long, mean)
factornames          <- c("ID", "taskOrder", "conditionOrder", "condition")
modelResults.summary[factornames] <- lapply(modelResults.summary[factornames], as.factor) 

fit.orderEffects <- aov(r_param ~ (conditionOrder * taskOrder * condition) + 
            Error(ID/condition) + conditionOrder * taskOrder,
            data = modelResults.summary)
summary(fit.orderEffects)
## 
## Error: ID
##                          Df Sum Sq Mean Sq F value Pr(>F)
## conditionOrder            1  0.012 0.01246   0.087  0.769
## taskOrder                 1  0.158 0.15830   1.103  0.296
## conditionOrder:taskOrder  1  0.002 0.00212   0.015  0.904
## Residuals                95 13.629 0.14346               
## 
## Error: ID:condition
##                                    Df Sum Sq Mean Sq F value Pr(>F)    
## condition                           1  6.069   6.069 102.088 <2e-16 ***
## conditionOrder:condition            1  0.114   0.114   1.912 0.1700    
## taskOrder:condition                 1  0.206   0.206   3.474 0.0654 .  
## conditionOrder:taskOrder:condition  1  0.157   0.157   2.634 0.1079    
## Residuals                          95  5.647   0.059                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.1.0.1 BF: Confirm no main effect of condition order on r_parameter

rParam.mainEffectCondAlone.bf       = lmBF(r_param ~ condition + ID,                  data=modelResults.long,  whichRandom="ID")
rParam.mainEffectsTaskOrder.bf      = lmBF(r_param ~ condition + taskOrder + ID,      data=modelResults.long,  whichRandom="ID")
rParam.mainEffectsEnvironOrder.bf   = lmBF(r_param ~ condition + conditionOrder + ID, data=modelResults.long,  whichRandom="ID")

# Compare the two models
bf.condOrder = rParam.mainEffectsEnvironOrder.bf / rParam.mainEffectCondAlone.bf
bf.condOrder
## Bayes factor analysis
## --------------
## [1] condition + conditionOrder + ID : 0.2296611 ±2.81%
## 
## Against denominator:
##   r_param ~ condition + ID 
## ---
## Bayes factor type: BFlinearModel, JZS
bf.taskOrder = rParam.mainEffectsTaskOrder.bf / rParam.mainEffectCondAlone.bf
bf.taskOrder
## Bayes factor analysis
## --------------
## [1] condition + taskOrder + ID : 0.3658266 ±2.09%
## 
## Against denominator:
##   r_param ~ condition + ID 
## ---
## Bayes factor type: BFlinearModel, JZS

1.1.1 Confirm no interaction of condition and conditionOrder on r_parameter

rParam.mainEffects.bf           = lmBF(r_param ~ condition + taskOrder + conditionOrder + ID,                            data=modelResults.long,  whichRandom="ID")
rParam.InteractionTaskOrder.bf  = lmBF(r_param ~ condition + taskOrder + conditionOrder + condition:taskOrder + ID,      data=modelResults.long,  whichRandom="ID")
rParam.InteractionCondOrder.bf  = lmBF(r_param ~ condition + taskOrder + conditionOrder + condition:conditionOrder + ID, data=modelResults.long,  whichRandom="ID")

# Compare the two models
bf.condXtaskOrder = rParam.InteractionTaskOrder.bf / rParam.mainEffects.bf
bf.condXtaskOrder
## Bayes factor analysis
## --------------
## [1] condition + taskOrder + conditionOrder + condition:taskOrder + ID : 0.8019386 ±5.58%
## 
## Against denominator:
##   r_param ~ condition + taskOrder + conditionOrder + ID 
## ---
## Bayes factor type: BFlinearModel, JZS
bf.condXCondOrder = rParam.InteractionCondOrder.bf / rParam.mainEffects.bf
bf.condXCondOrder
## Bayes factor analysis
## --------------
## [1] condition + taskOrder + conditionOrder + condition:conditionOrder + ID : 0.2791764 ±4.63%
## 
## Against denominator:
##   r_param ~ condition + taskOrder + conditionOrder + ID 
## ---
## Bayes factor type: BFlinearModel, JZS
# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

2 Congruency index

Compare indeces between conditions, partially to show that the manipulation worked ### Population

sd(congruencyIndex$population)
## [1] 0.09665959
t.test(congruencyIndex$population, mu=0.5, var.equal=TRUE)
## 
##  One Sample t-test
## 
## data:  congruencyIndex$population
## t = 30.569, df = 98, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0.5
## 95 percent confidence interval:
##  0.7776913 0.8162481
## sample estimates:
## mean of x 
## 0.7969697
cohensD(congruencyIndex$population, mu=0.5)             # no method = "paired" in this case: one sample tests!
## [1] 3.072325

2.0.1 Distance

sd(congruencyIndex$distance)
## [1] 0.1332554
t.test(congruencyIndex$distance , mu=0.5, var.equal=TRUE)
## 
##  One Sample t-test
## 
## data:  congruencyIndex$distance
## t = 14.61, df = 98, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0.5
## 95 percent confidence interval:
##  0.6690937 0.7222483
## sample estimates:
## mean of x 
##  0.695671
cohensD(congruencyIndex$distance, mu=0.5) 
## [1] 1.468391
# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

3 Lost Focus

3.1 Identify distracted subjects during recognition task

distractedSubjects_Recognition <- sum(resultsRH.perSubj$totalLostFocusRecognized)
distractedSubjects_Recognition
## [1] 6

3.2 Identify distracted subjects during comparative judgement task

distractedSubjects <- aggregate(lostFocus ~ subject, resultsRH.perTrial.all, sum)
taskAccuracy       <- aggregate(correct ~ subject, resultsRH.perTrial.all, mean)

3.3 Check overlap between distracted subjects in recognition and comparative judgement task

distractedSubjects$lostFocus[(resultsRH.perSubj$totalLostFocusRecognized>0)]
## [1]  1 17  4  2 11 10

3.3.1 Check whether mean task accuracy was higher for more distracted subjects

cor.test(distractedSubjects$lostFocus, taskAccuracy$correct) 
## 
##  Pearson's product-moment correlation
## 
## data:  distractedSubjects$lostFocus and taskAccuracy$correct
## t = 0.036937, df = 97, p-value = 0.9706
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1938050  0.2010133
## sample estimates:
##         cor 
## 0.003750317
t.test(taskAccuracy$correct[distractedSubjects$lostFocus >= 1], taskAccuracy$correct[distractedSubjects$lostFocus < 1])
## 
##  Welch Two Sample t-test
## 
## data:  taskAccuracy$correct[distractedSubjects$lostFocus >= 1] and taskAccuracy$correct[distractedSubjects$lostFocus < 1]
## t = -1.468, df = 96.825, p-value = 0.1454
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.043543893  0.006517654
## sample estimates:
## mean of x mean of y 
## 0.6314869 0.6500000
ttestBF(x = taskAccuracy$correct[distractedSubjects$lostFocus >= 1], y = taskAccuracy$correct[distractedSubjects$lostFocus < 1], paired=FALSE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.549162 ±0%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS
distracted <- data.frame(lostFocusTimes = distractedSubjects$lostFocus,
                         meanAccuracy = taskAccuracy$correct)
lmBF(meanAccuracy ~ lostFocusTimes, data = distracted)
## Bayes factor analysis
## --------------
## [1] lostFocusTimes : 0.2118705 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

3.4 So we decide to remove all subjects that left the browser window at least once

toRemove_index     <- which(distractedSubjects$lostFocus >= 1)
toRemove_subject   <- distractedSubjects$subject[distractedSubjects$lostFocus >= 1]
length(toRemove_subject)
## [1] 49
# Clean up data
congruencyIndex         <- congruencyIndex[-toRemove_index, ] 
modelResults            <- modelResults[-toRemove_index, ]
modelResults.long       <- modelResults.long[-c(modelResults.long$ID %in% toRemove_index), ]
resultsRH.perSubj       <- resultsRH.perSubj[-toRemove_index, ]
resultsRH.perTrial.all  <- resultsRH.perTrial.all[-c(which(resultsRH.perTrial.all$subject %in% toRemove_subject)), ] 

# Run all the analyses again
# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::


library(lsr)        # compute effect sizes (cohensD) in ttests
library(R.matlab)
library(plyr)       # to summarize data and required before dplyr
library(dplyr)      # for renaming columns 
library(reshape)    # melt and cast functions
source(paste(fullPath,"/functions/as.data.frame.list.R", sep=""))
library(BayesFactor)

# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

4 Criterion knowledge

# Melt for ggplot2
criterionKnowledge_toPlot <- resultsRH.perSubj[,c("rankingRhoRecogOnly_pop", "rankingRhoRecogOnly_dist",
                                                                          "rankingRho_pop",              "rankingRho_dist")]
criterionKnowledge_toPlot <- melt(criterionKnowledge_toPlot)
## Using  as id variables
# Just some renaming for convenience 
# assign condition ID
criterionKnowledge_toPlot$condition[criterionKnowledge_toPlot$variable == "rankingRhoRecogOnly_pop"]    <- "pop"
criterionKnowledge_toPlot$condition[criterionKnowledge_toPlot$variable == "rankingRho_pop"]                 <- "pop"
criterionKnowledge_toPlot$condition[criterionKnowledge_toPlot$variable == "rankingRhoRecogOnly_dist"] <- "dist"
criterionKnowledge_toPlot$condition[criterionKnowledge_toPlot$variable == "rankingRho_dist"]                <- "dist"
criterionKnowledge_toPlot$condition <- as.factor(criterionKnowledge_toPlot$condition)
# assign "Recognized only" ID
criterionKnowledge_toPlot$Ronly[criterionKnowledge_toPlot$variable == "rankingRhoRecogOnly_pop"]        <- "Ronly"
criterionKnowledge_toPlot$Ronly[criterionKnowledge_toPlot$variable == "rankingRho_pop"]                     <- "allcities"
criterionKnowledge_toPlot$Ronly[criterionKnowledge_toPlot$variable == "rankingRhoRecogOnly_dist"]     <- "Ronly"
criterionKnowledge_toPlot$Ronly[criterionKnowledge_toPlot$variable == "rankingRho_dist"]                    <- "allcities"
criterionKnowledge_toPlot$Ronly <- as.factor(criterionKnowledge_toPlot$Ronly)

library(ggplot2)
ggplot(criterionKnowledge_toPlot, aes(x = condition, y = value, color = condition)) +
        geom_jitter() +
        facet_grid(. ~ Ronly) +
        ylab("criterionKnowledge_toPlot")

rhoCritKnow.mean <- aggregate(value ~ condition + Ronly, data = criterionKnowledge_toPlot, mean)
rhoCritKnow.mean
##   condition     Ronly     value
## 1      dist allcities 0.1924999
## 2       pop allcities 0.5110001
## 3      dist     Ronly 0.3925503
## 4       pop     Ronly 0.3705307
rhoCritKnow.sd   <- meanVals <- aggregate(value ~ condition + Ronly, data = criterionKnowledge_toPlot, sd)
rhoCritKnow.sd
##   condition     Ronly     value
## 1      dist allcities 0.3465933
## 2       pop allcities 0.1845386
## 3      dist     Ronly 0.4543289
## 4       pop     Ronly 0.3684305

4.1 Stats on criterion Knowledge

4.1.1 Permutation test for the recognized-only cities

library(exactRankTests)
perm.test(resultsRH.perSubj$rankingRhoRecogOnly_pop, resultsRH.perSubj$rankingRhoRecogOnly_dist, paired=TRUE, alternative="two.sided")
## 
##  1-sample Permutation Test (scores mapped into 1:m using rounded
##  scores)
## 
## data:  resultsRH.perSubj$rankingRhoRecogOnly_pop and resultsRH.perSubj$rankingRhoRecogOnly_dist
## T = 231, p-value = 0.09777
## alternative hypothesis: true mu is not equal to 0

4.1.2 BF

rhoCritKnow.bf = ttestBF(x = resultsRH.perSubj$rankingRhoRecogOnly_pop,
                         y = resultsRH.perSubj$rankingRhoRecogOnly_dist, paired=TRUE)
rhoCritKnow.bf
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.1602746 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS

4.1.3 Permutation test for the all cities

perm.test(resultsRH.perSubj$rankingRho_pop, resultsRH.perSubj$rankingRho_dist, paired=TRUE, alternative="two.sided")
## 
##  1-sample Permutation Test (scores mapped into 1:m using rounded
##  scores)
## 
## data:  resultsRH.perSubj$rankingRho_pop and resultsRH.perSubj$rankingRho_dist
## T = 741, p-value = 6.592e-08
## alternative hypothesis: true mu is not equal to 0

4.1.4 BF

rhoCritKnow.bf = ttestBF(x = resultsRH.perSubj$rankingRho_pop,
             y = resultsRH.perSubj$rankingRho_dist, paired=TRUE)
rhoCritKnow.bf
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 86687.49 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

5 Recognition Rate and Applicability of the RH

recognized.percent.mean <- mean(resultsRH.perSubj$totalRecognized/15*100)
recognized.percent.mean
## [1] 52.66667
recognized.percent.sd   <- sd(resultsRH.perSubj$totalRecognized/15*100)
recognized.percent.sd
## [1] 11.60577

5.1 Proportion of RU trials in the Comparative Judgment Task

resultsRH.perTrial.all$subject <- as.factor(resultsRH.perTrial.all$subject)
trialProportion <- resultsRH.perTrial.all[resultsRH.perTrial.all$condition==1,] %>% #only one condition is necessary. It's the same for both
    group_by(subject) %>% 
    summarise(
      UUtrials  = length(which(trialType==0)),
      RUtrials  = length(which(trialType==1)),
      RRtrials  = length(which(trialType==2))
      )
trialProportion <- as.data.frame(trialProportion)

# Check output: should all be 70 (sum all columns except subject)
all(rowSums(trialProportion[,!names(trialProportion) %in% c("subject")]) == 70)
## [1] TRUE
# Check output: should all be 100 (sum all columns except subject)
trialProportion.percentage <- trialProportion[,!names(trialProportion) %in% c("subject")]/70*100
all(round(rowSums(trialProportion.percentage[,!names(trialProportion.percentage) %in% c("subject")])) == 100)
## [1] TRUE
mean(trialProportion.percentage$RUtrials)
## [1] 50.74286
sd(trialProportion.percentage$RUtrials)
## [1] 7.02373

5.1.1 Plot all subjs’ trial proportions

trialProportion.sorted <- trialProportion[order(trialProportion$UUtrials),]
trialProportion.sorted$subj.sorted <- c(1:dim(trialProportion.sorted)[1])
trialProportion.sorted <- reshape2:::melt.data.frame(trialProportion.sorted,  id.vars=c("subj.sorted"),    #explicitly call melt from reshape2 otherwise coumns might not get properly named
                                                        measure.vars = c("UUtrials", "RUtrials", "RRtrials"),
                                                        variable.name = "trialType",
                                                        value.name = "trialCount")
ggplot(trialProportion.sorted, aes(x = subj.sorted, y = trialCount, fill=trialType)) +
    geom_bar(stat='identity')

# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

6 (Performance in the) Comparative Judgment Task

6.1 Population vs distance

# Group by subject
proportionCorrect.perSubj <- resultsRH.perTrial.all %>% 
    group_by(subject, condition) %>% 
    summarise(
      correct  = mean(correct, na.rm = TRUE)
      )

proportionCorrect.perSubj <- as.data.frame(proportionCorrect.perSubj)
proportionCorrect.perSubj$condition[proportionCorrect.perSubj$condition == 1] <- "pop"
proportionCorrect.perSubj$condition[proportionCorrect.perSubj$condition == 2] <- "dist"

proportionCorrect.perSubj.wide <- reshape(proportionCorrect.perSubj, timevar=c("condition"), 
                  idvar=c("subject"), dir="wide")

6.1.1 ttest

t.test(proportionCorrect.perSubj.wide$correct.pop, proportionCorrect.perSubj.wide$correct.dist, paired = TRUE)
## 
##  Paired t-test
## 
## data:  proportionCorrect.perSubj.wide$correct.pop and proportionCorrect.perSubj.wide$correct.dist
## t = 3.739, df = 49, p-value = 0.0004838
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.02748825 0.09136890
## sample estimates:
## mean of the differences 
##              0.05942857
cohensD(proportionCorrect.perSubj.wide$correct.pop, proportionCorrect.perSubj.wide$correct.dist, method="paired")
## [1] 0.5287802

6.1.2 BF

proportionCorrect.bf = ttestBF(x = proportionCorrect.perSubj.wide$correct.pop,
                               y = proportionCorrect.perSubj.wide$correct.dist, paired=TRUE)
proportionCorrect.bf
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 54.70244 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS

6.2 Compare population vs. distance

6.2.1 alphas

t.test(modelResults$alpha.1, modelResults$alpha.2, paired=TRUE, var.equal=TRUE)  #Because of renaming, 1 is population and 2 is distance
## 
##  Paired t-test
## 
## data:  modelResults$alpha.1 and modelResults$alpha.2
## t = 16.28, df = 49, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2860461 0.3666075
## sample estimates:
## mean of the differences 
##               0.3263268
cohensD(modelResults$alpha.1, modelResults$alpha.2, method = "paired") 
## [1] 2.302367
ttestBF(x = modelResults$alpha.1, y = modelResults$alpha.2, paired=TRUE)
## t is large; approximation invoked.
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 2.085345e+18 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS

6.2.2 r parameters

t.test(modelResults$r_param.1, modelResults$r_param.2, paired=TRUE, var.equal=TRUE)
## 
##  Paired t-test
## 
## data:  modelResults$r_param.1 and modelResults$r_param.2
## t = 5.702, df = 49, p-value = 6.714e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2030770 0.4241243
## sample estimates:
## mean of the differences 
##               0.3136006
cohensD(modelResults$r_param.1, modelResults$r_param.2, method = "paired") 
## [1] 0.8063822
ttestBF(x = modelResults$r_param.1, y = modelResults$r_param.2, paired=TRUE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 23786.09 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS

6.2.3 b parameters

t.test(modelResults$b_param.1, modelResults$b_param.2, paired=TRUE, var.equal=TRUE)  #Because of renaming, 1 is population and 2 is distance
## 
##  Paired t-test
## 
## data:  modelResults$b_param.1 and modelResults$b_param.2
## t = -3.5124, df = 49, p-value = 0.0009648
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.07321773 -0.01992649
## sample estimates:
## mean of the differences 
##             -0.04657211
cohensD(modelResults$b_param.1, modelResults$b_param.2, method = "paired") 
## [1] 0.4967289
# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

7 Formal Modeling: Adaptive Use of the RH

This section is not included here

# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

8 Metacognitive Sensitivity

Do a glm approach, as described my Siedlecka et al (2016, Front Psychol)

library(lme4)
library(lmerTest)

# Prepare the dataset
# Separate between RU+ and RU-
resultsRH.perTrial.all$trialType <- rowSums(cbind(resultsRH.perTrial.all$trialType,0.5 *resultsRH.perTrial.all$RU.Rchosen), na.rm=TRUE)
# Check output
unique(resultsRH.perTrial.all$trialType)
## [1] 0.0 1.5 2.0 0.5
# Set factors
resultsRH.perTrial.all$trialType <- as.factor(resultsRH.perTrial.all$trialType)
resultsRH.perTrial.all$condition <- as.factor(resultsRH.perTrial.all$condition)

# Add the sorted cityPair
# But be careful with this, because city pair will be (partially) confounded with condition
twoCities <- resultsRH.perTrial.all[,c("city1", "city2")]
twoCities <- as.data.frame(t(apply(twoCities, 1, sort))) #Sort the two cities in alphabetical order (so no matter how they were shown on the screen)
resultsRH.perTrial.all$cityPair <- as.factor(paste(twoCities$V1,twoCities$V2,sep=""))


# Get only RU+ trials 
myDataForglm.RU                   <- subset(resultsRH.perTrial.all, trialType %in% c("1.5")) 
myDataForglm.RU                   <- as.data.frame(myDataForglm.RU)
myDataForglm.RU$trialType         <- factor(myDataForglm.RU$trialType)         
myDataForglm.RU$choiceConf.scaled <- scale(myDataForglm.RU$choiceConf)

8.1 Do the models on RU+ trials only

glmres.RU.conf.full = glmer(correct ~ choiceConf.scaled * condition + (1 + choiceConf.scaled|subject)  + (1|cityPair), data = myDataForglm.RU, family = binomial)
glmres.RU.conf.sum  = glmer(correct ~ choiceConf.scaled + condition + (1 + choiceConf.scaled|subject)  + (1|cityPair), data = myDataForglm.RU, family = binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0289017 (tol =
## 0.001, component 1)
# The full model fails to converge but the subsequent tests show that the obtained result is actually reasonable
source(system.file("utils", "allFit.R", package="lme4"))
## Loading required package: optimx
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'optimx'
## Loading required package: dfoptim
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'dfoptim'
fm1.all <- allFit(glmres.RU.conf.full)
## bobyqa : [OK]
## Nelder_Mead :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00197747 (tol =
## 0.001, component 1)
## [OK]
## nlminbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
## [OK]
## nloptwrap.NLOPT_LN_BOBYQA :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
## [OK]
ss <- summary(fm1.all)
ss$fixef
##                               (Intercept) choiceConf.scaled condition2
## bobyqa                           2.300971        0.03740569  -1.915330
## Nelder_Mead                      2.300976        0.03740393  -1.915336
## nlminbw                          2.300963        0.03740474  -1.915333
## nloptwrap.NLOPT_LN_NELDERMEAD    2.251041        0.04603965  -1.895663
## nloptwrap.NLOPT_LN_BOBYQA        2.251041        0.04603965  -1.895663
##                               choiceConf.scaled:condition2
## bobyqa                                           0.1588917
## Nelder_Mead                                      0.1588941
## nlminbw                                          0.1588909
## nloptwrap.NLOPT_LN_NELDERMEAD                    0.1508891
## nloptwrap.NLOPT_LN_BOBYQA                        0.1508891

8.1.1 χ^2

interactionModelComp <- anova(glmres.RU.conf.full,glmres.RU.conf.sum) 
interactionModelComp
## Data: myDataForglm.RU
## Models:
## glmres.RU.conf.sum: correct ~ choiceConf.scaled + condition + (1 + choiceConf.scaled | 
## glmres.RU.conf.sum:     subject) + (1 | cityPair)
## glmres.RU.conf.full: correct ~ choiceConf.scaled * condition + (1 + choiceConf.scaled | 
## glmres.RU.conf.full:     subject) + (1 | cityPair)
##                     Df  AIC    BIC logLik deviance  Chisq Chi Df
## glmres.RU.conf.sum   7 2398 2439.5  -1192     2384              
## glmres.RU.conf.full  8 2398 2445.3  -1191     2382 2.0303      1
##                     Pr(>Chisq)
## glmres.RU.conf.sum            
## glmres.RU.conf.full     0.1542

8.1.2 BF10

interactionModelComp.BF10 = 1/exp((interactionModelComp[3][[1]][2] - interactionModelComp[3][[1]][1])/2) #BF10
interactionModelComp.BF10
## [1] 0.05265585
# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

9 Individual confidence ratings

meanConfPerSubj <- ddply(resultsRH.perTrial.all, c("trialType", "condition", "subject"), summarise,
                          meanConf = mean(choiceConf)
                        )
meanConfPerSubj.RUplus      <- subset(meanConfPerSubj, trialType %in% c(1.5))
meanConfPerSubj.RUplus.pop  <- subset(meanConfPerSubj.RUplus, condition %in% c(1)) 
meanConfPerSubj.RUplus.dist <- subset(meanConfPerSubj.RUplus, condition %in% c(2)) 

9.1 α and mean confidence in RU+

forBF.alphaVsConf <- data.frame(alpha.pop            = modelResults$alpha.1,
                                alpha.dist           = modelResults$alpha.2,
                                meanConf.RUplus.pop  = meanConfPerSubj.RUplus.pop$meanConf,
                                meanConf.RUplus.dist = meanConfPerSubj.RUplus.dist$meanConf,
                                ID                   = modelResults$ID)

# ### population
cor.test(forBF.alphaVsConf$alpha.pop, forBF.alphaVsConf$meanConf.RUplus.pop)
## 
##  Pearson's product-moment correlation
## 
## data:  forBF.alphaVsConf$alpha.pop and forBF.alphaVsConf$meanConf.RUplus.pop
## t = -0.32841, df = 48, p-value = 0.744
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3214604  0.2340835
## sample estimates:
##         cor 
## -0.04734934
lmBF(meanConf.RUplus.pop ~ alpha.pop + ID, data = forBF.alphaVsConf)
## Bayes factor analysis
## --------------
## [1] alpha.pop + ID : 0.3692062 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
# ### distance
cor.test(forBF.alphaVsConf$alpha.dist, forBF.alphaVsConf$meanConf.RUplus.dist)
## 
##  Pearson's product-moment correlation
## 
## data:  forBF.alphaVsConf$alpha.dist and forBF.alphaVsConf$meanConf.RUplus.dist
## t = -1.44, df = 48, p-value = 0.1564
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.45601437  0.07934549
## sample estimates:
##        cor 
## -0.2034966
lmBF(meanConf.RUplus.dist ~ alpha.dist + ID, data = forBF.alphaVsConf)
## Bayes factor analysis
## --------------
## [1] alpha.dist + ID : 0.5139567 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

9.2 Correlation between Δconfidence and Δα: (∆pop-dist)

RU.Delta <- data.frame(meanConf_RH = meanConfPerSubj.RUplus.pop$meanConf - meanConfPerSubj.RUplus.dist$meanConf,
                       alpha       = modelResults$alpha.1                - modelResults$alpha.2,
                       ID          = modelResults$ID)

cor.test(RU.Delta$alpha, RU.Delta$meanConf_RH)
## 
##  Pearson's product-moment correlation
## 
## data:  RU.Delta$alpha and RU.Delta$meanConf_RH
## t = -0.24579, df = 48, p-value = 0.8069
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3107352  0.2453146
## sample estimates:
##         cor 
## -0.03545403
# plot(RU.Delta$alpha, RU.Delta$meanConf_RH)
lmBF(meanConf_RH ~ alpha + ID, data = RU.Delta)
## Bayes factor analysis
## --------------
## [1] alpha + ID : 0.1173714 ±0.01%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

10 Mean confidence ratings

10.1 ANOVA with all trial types: RU+, RU- (4x2)

3x2 anova on mean confidence with all trial types (RR, RUplus, RUminus, UU) * condition

#  Get means
myDataForglm.RUp.RUm.RR.UU                   <- subset(resultsRH.perTrial.all, trialType %in% c("0", "0.5", "1.5", "2")) 
myDataForglm.RUp.RUm.RR.UU                   <- as.data.frame(myDataForglm.RUp.RUm.RR.UU)
myDataForglm.RUp.RUm.RR.UU$trialType         <- factor(myDataForglm.RUp.RUm.RR.UU$trialType)        
myDataForglm.RUp.RUm.RR.UU$choiceConf.scaled <- scale(myDataForglm.RUp.RUm.RR.UU$choiceConf)
myDataForglm.RUp.RUm.RR.UU.mean              <- aggregate(choiceConf ~ subject + condition + trialType, myDataForglm.RUp.RUm.RR.UU, mean)


#Now look for subjects that did not have any RU- or no UU trials and remove them
yesUU.pop   <- unique(myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$trialType == 0   & myDataForglm.RUp.RUm.RR.UU.mean$condition == 1])
yesUU.dist  <- unique(myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$trialType == 0   & myDataForglm.RUp.RUm.RR.UU.mean$condition == 2])
yesRUm.pop  <- unique(myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$trialType == 0.5 & myDataForglm.RUp.RUm.RR.UU.mean$condition == 1])
yesRUm.dist <- unique(myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$trialType == 0.5 & myDataForglm.RUp.RUm.RR.UU.mean$condition == 2])
yesRUp.pop  <- unique(myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$trialType == 1.5 & myDataForglm.RUp.RUm.RR.UU.mean$condition == 1])
yesRUp.dist <- unique(myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$trialType == 1.5 & myDataForglm.RUp.RUm.RR.UU.mean$condition == 2])
yesRR       <- unique(myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$trialType == 2])
noUU.pop    <- as.numeric(setdiff(yesRR,yesUU.pop))
noUU.dist   <- as.numeric(setdiff(yesRR,yesUU.dist))
noRUm.pop   <- as.numeric(setdiff(yesRR,yesRUm.pop))
noRUm.dist  <- as.numeric(setdiff(yesRR,yesRUm.dist))
noRUp.pop   <- as.numeric(setdiff(yesRR,yesRUp.pop))
noRUp.dist  <- as.numeric(setdiff(yesRR,yesRUp.dist))
missings <- unique(c(noUU.pop, noUU.dist, noRUm.pop, noRUm.dist, noRUp.pop, noRUp.dist))
myDataForglm.RUp.RUm.RR.UU.mean$subject[myDataForglm.RUp.RUm.RR.UU.mean$subject %in% missings]  <- NA
myDataForglm.RUp.RUm.RR.UU.mean <- myDataForglm.RUp.RUm.RR.UU.mean[complete.cases(myDataForglm.RUp.RUm.RR.UU.mean),]

factornames <- c("subject", "condition", "trialType")
myDataForglm.RUp.RUm.RR.UU.mean[factornames] <- lapply(myDataForglm.RUp.RUm.RR.UU.mean[factornames], as.factor) 

fit.meanConfs.RUp.RUm.RR.UU <- aov(
            choiceConf ~ (condition * trialType) + 
            Error(subject/(condition * trialType)),
            data = myDataForglm.RUp.RUm.RR.UU.mean)
summary(fit.meanConfs.RUp.RUm.RR.UU)
## 
## Error: subject
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 33  81282    2463               
## 
## Error: subject:condition
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## condition  1    985   984.9   7.843 0.00847 **
## Residuals 33   4144   125.6                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: subject:trialType
##           Df Sum Sq Mean Sq F value Pr(>F)    
## trialType  3  41723   13908   61.72 <2e-16 ***
## Residuals 99  22307     225                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: subject:condition:trialType
##                     Df Sum Sq Mean Sq F value Pr(>F)
## condition:trialType  3    310  103.26   1.272  0.288
## Residuals           99   8038   81.19

10.1.1 ANOVA effect sizes

library(DescTools)
EtaSq(fit.meanConfs.RUp.RUm.RR.UU, type=1, anova=TRUE)
##                          eta.sq eta.sq.part  eta.sq.gen         SS df
## condition           0.006202387   0.1920211 0.008435253   984.8668  1
## trialType           0.262755633   0.6516170 0.264915550 41722.5312  3
## condition:trialType 0.001950842   0.0371081 0.002668577   309.7709  3
##                             MS       SSE dfE         F            p
## condition             984.8668  4144.085  33  7.842648 8.465095e-03
## trialType           13907.5104 22306.697  99 61.723325 1.393558e-22
## condition:trialType   103.2570  8038.025  99  1.271760 2.882939e-01

10.2 ANOVA with all trial types: RU regardless of +/- (3x2)

3x2 anova on mean confidence with all trial types (RR, RU, UU) * condition

#  Get means
myDataForglm.RU.RR.UU                   <- subset(resultsRH.perTrial.all, trialType %in% c("0", "1.5", "2")) 
myDataForglm.RU.RR.UU                   <- as.data.frame(myDataForglm.RU.RR.UU)
myDataForglm.RU.RR.UU$trialType         <- factor(myDataForglm.RU.RR.UU$trialType)        
myDataForglm.RU.RR.UU$choiceConf.scaled <- scale(myDataForglm.RU.RR.UU$choiceConf)

10.3 ANOVA only with RR, RU+

myDataForglm.RU.RR.UU.mean <- aggregate(choiceConf ~ subject + condition + trialType, myDataForglm.RU.RR.UU, mean)

#controlcount <- aggregate(choiceConf ~ subject + condition + trialType, myDataForglm.RU.RR.UU, length)
#Subject 45 did not have any UU trials in one of the conditions. So I remove it 

myDataForglm.RU.RR.UU.mean$choiceConf[myDataForglm.RU.RR.UU.mean$subject == 45] <- NA 
myDataForglm.RU.RR.UU.mean <- myDataForglm.RU.RR.UU.mean[complete.cases(myDataForglm.RU.RR.UU.mean),]
factornames <- c("subject", "condition", "trialType")
myDataForglm.RU.RR.UU.mean[factornames] <- lapply(myDataForglm.RU.RR.UU.mean[factornames], as.factor) 

fit.meanConfs.RU.RR.UU <- aov(
            choiceConf ~ (condition * trialType) + 
            Error(subject/(condition * trialType)),
            data = myDataForglm.RU.RR.UU.mean)
summary(fit.meanConfs.RU.RR.UU)
## 
## Error: subject
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 48  76649    1597               
## 
## Error: subject:condition
##           Df Sum Sq Mean Sq F value  Pr(>F)   
## condition  1   1798  1797.6    9.91 0.00283 **
## Residuals 48   8707   181.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: subject:trialType
##           Df Sum Sq Mean Sq F value Pr(>F)    
## trialType  2  54565   27283   116.1 <2e-16 ***
## Residuals 96  22563     235                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: subject:condition:trialType
##                     Df Sum Sq Mean Sq F value  Pr(>F)   
## condition:trialType  2    886   443.1   5.877 0.00391 **
## Residuals           96   7238    75.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

10.3.1 ANOVA effect sizes

library(DescTools)
EtaSq(fit.meanConfs.RU.RR.UU, type=1, anova=TRUE)
##                          eta.sq eta.sq.part  eta.sq.gen         SS df
## condition           0.010426644   0.1711237 0.015370168  1797.6291  1
## trialType           0.316490997   0.7074609 0.321495635 54565.3424  2
## condition:trialType 0.005140092   0.1090753 0.007636647   886.1891  2
##                             MS       SSE dfE          F            p
## condition            1797.6291  8707.220  48   9.909729 2.825330e-03
## trialType           27282.6712 22563.080  96 116.080624 2.381382e-26
## condition:trialType   443.0946  7238.373  96   5.876607 3.911721e-03

2x2 anova on mean confidence with trial type (RU+ vs RR) * condition

#  Get means
# myDataForglm.RU.RR                   <- subset(resultsRH.perTrial.all, trialType %in% c("1.5", "2")) 
# myDataForglm.RU.RR                   <- as.data.frame(myDataForglm.RU.RR)
# myDataForglm.RU.RR$trialType         <- factor(myDataForglm.RU.RR$trialType)        
# myDataForglm.RU.RR$choiceConf.scaled <- scale(myDataForglm.RU.RR$choiceConf)

# #' ## Anova
# myDataForglm.RU.RR.mean <- aggregate(choiceConf ~ subject + condition + trialType, myDataForglm.RU.RR, mean)
# factornames <- c("subject", "condition", "trialType")
# myDataForglm.RU.RR.mean[factornames] <- lapply(myDataForglm.RU.RR.mean[factornames], as.factor) 

# fit.meanConfs <- aov(
#             choiceConf ~ (condition * trialType) + 
#             Error(subject/(condition * trialType)),
#             data = myDataForglm.RU.RR.mean)
# summary(fit.meanConfs)

# #' ### ANOVA effect sizes
# library(DescTools)
# EtaSq(fit.meanConfs, type=1, anova=TRUE)

# # bayesian anova
# # anovaBF(choiceConf ~ condition * trialType + subject, data = myDataForglm.RU.RR.mean,  whichRandom="subject")

# #' ### BF for interaction
# choiceConf.sum.bf  = lmBF(choiceConf ~ condition + trialType + subject,                       data=myDataForglm.RU.RR.mean,  whichRandom="subject")
# choiceConf.full.bf = lmBF(choiceConf ~ condition + trialType + condition:trialType + subject, data=myDataForglm.RU.RR.mean,  whichRandom="subject")

# # Compare the two models
# bf.condXtrialType = choiceConf.full.bf / choiceConf.sum.bf
# bf.condXtrialType


# #' ## Follow-up t-tests
# # Get each condition (probably could have been done with aggregate)
# RU.pop.means  <- myDataForglm.RU.RR.mean[myDataForglm.RU.RR.mean$trialType=="1.5" & myDataForglm.RU.RR.mean$condition=="1",]
# RR.pop.means  <- myDataForglm.RU.RR.mean[myDataForglm.RU.RR.mean$trialType=="2"   & myDataForglm.RU.RR.mean$condition=="1",]
# RU.dist.means <- myDataForglm.RU.RR.mean[myDataForglm.RU.RR.mean$trialType=="1.5" & myDataForglm.RU.RR.mean$condition=="2",]
# RR.dist.means <- myDataForglm.RU.RR.mean[myDataForglm.RU.RR.mean$trialType=="2"   & myDataForglm.RU.RR.mean$condition=="2",]




RR.pop.means     <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="2"   & myDataForglm.RUp.RUm.RR.UU.mean$condition=="1",]
RR.dist.means    <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="2"   & myDataForglm.RUp.RUm.RR.UU.mean$condition=="2",]
RUplu.pop.means  <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="1.5" & myDataForglm.RUp.RUm.RR.UU.mean$condition=="1",]
RUplu.dist.means <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="1.5" & myDataForglm.RUp.RUm.RR.UU.mean$condition=="2",]
RUmin.pop.means  <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="0.5" & myDataForglm.RUp.RUm.RR.UU.mean$condition=="1",]
RUmin.dist.means <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="0.5" & myDataForglm.RUp.RUm.RR.UU.mean$condition=="2",]
UU.pop.means     <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="0"   & myDataForglm.RUp.RUm.RR.UU.mean$condition=="1",]
UU.dist.means    <- myDataForglm.RUp.RUm.RR.UU.mean[myDataForglm.RUp.RUm.RR.UU.mean$trialType=="0"   & myDataForglm.RUp.RUm.RR.UU.mean$condition=="2",]

10.3.2 RR trials

t.test(RR.pop.means$choiceConf, RR.dist.means$choiceConf, paired = TRUE)
## 
##  Paired t-test
## 
## data:  RR.pop.means$choiceConf and RR.dist.means$choiceConf
## t = 0.49202, df = 33, p-value = 0.626
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.484887  5.708105
## sample estimates:
## mean of the differences 
##                1.111609
cohensD(RR.pop.means$choiceConf, RR.dist.means$choiceConf, method="paired")
## [1] 0.08438141
ttestBF(x = RR.pop.means$choiceConf, y = RR.dist.means$choiceConf, paired=TRUE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.2056337 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
mean(RR.pop.means$choiceConf - RR.dist.means$choiceConf)
## [1] 1.111609
sd(RR.pop.means$choiceConf - RR.dist.means$choiceConf)
## [1] 13.17363

10.3.3 RU+ trials

t.test(RUplu.pop.means$choiceConf, RUplu.dist.means$choiceConf, paired = TRUE)
## 
##  Paired t-test
## 
## data:  RUplu.pop.means$choiceConf and RUplu.dist.means$choiceConf
## t = 3.0543, df = 33, p-value = 0.004439
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   2.268887 11.321635
## sample estimates:
## mean of the differences 
##                6.795261
cohensD(RUplu.pop.means$choiceConf, RUplu.dist.means$choiceConf, method="paired")
## [1] 0.5238141
ttestBF(x = RUplu.pop.means$choiceConf, y = RUplu.dist.means$choiceConf, paired=TRUE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 8.665426 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
mean(RUplu.pop.means$choiceConf - RUplu.dist.means$choiceConf)
## [1] 6.795261
sd(RUplu.pop.means$choiceConf - RUplu.dist.means$choiceConf)
## [1] 12.97266

10.3.4 RU- trials

t.test(RUmin.pop.means$choiceConf, RUmin.dist.means$choiceConf, paired = TRUE)
## 
##  Paired t-test
## 
## data:  RUmin.pop.means$choiceConf and RUmin.dist.means$choiceConf
## t = 1.012, df = 33, p-value = 0.3189
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.690254  8.014912
## sample estimates:
## mean of the differences 
##                2.662329
cohensD(RUmin.pop.means$choiceConf, RUmin.dist.means$choiceConf, method="paired")
## [1] 0.1735481
ttestBF(x = RUmin.pop.means$choiceConf, y = RUmin.dist.means$choiceConf, paired=TRUE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.2943693 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
mean(RUmin.pop.means$choiceConf - RUmin.dist.means$choiceConf)
## [1] 2.662329
sd(RUmin.pop.means$choiceConf - RUmin.dist.means$choiceConf)
## [1] 15.34058

10.3.5 UU trials

t.test(UU.pop.means$choiceConf, UU.dist.means$choiceConf, paired = TRUE)
## 
##  Paired t-test
## 
## data:  UU.pop.means$choiceConf and UU.dist.means$choiceConf
## t = 2.1376, df = 33, p-value = 0.04005
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2243866 9.0827961
## sample estimates:
## mean of the differences 
##                4.653591
cohensD(UU.pop.means$choiceConf, UU.dist.means$choiceConf, method="paired")
## [1] 0.3665928
ttestBF(x = UU.pop.means$choiceConf, y = UU.dist.means$choiceConf, paired=TRUE)
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 1.362484 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS
mean(UU.pop.means$choiceConf - UU.dist.means$choiceConf)
## [1] 4.653591
sd(UU.pop.means$choiceConf - UU.dist.means$choiceConf)
## [1] 12.69417

10.3.6 Difference (between trial types) of differences (between conditions)

# t.test((RU.pop.means$choiceConf - RU.dist.means$choiceConf), (RR.pop.means$choiceConf - RR.dist.means$choiceConf), paired = TRUE)
# cohensD((RU.pop.means$choiceConf - RU.dist.means$choiceConf), (RR.pop.means$choiceConf - RR.dist.means$choiceConf), method="paired")
# ttestBF(x = (RU.pop.means$choiceConf - RU.dist.means$choiceConf), y = (RR.pop.means$choiceConf - RR.dist.means$choiceConf), paired=TRUE)
# mean((RU.pop.means$choiceConf - RU.dist.means$choiceConf) - (RR.pop.means$choiceConf - RR.dist.means$choiceConf))
# sd((RU.pop.means$choiceConf - RU.dist.means$choiceConf) - (RR.pop.means$choiceConf - RR.dist.means$choiceConf))



# # ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

11 Table 2

11.0.1 Proportion correct

mean(proportionCorrect.perSubj.wide$correct.pop)
## [1] 0.6797143
mean(proportionCorrect.perSubj.wide$correct.dist)
## [1] 0.6202857
sd(proportionCorrect.perSubj.wide$correct.pop)
## [1] 0.06854774
sd(proportionCorrect.perSubj.wide$correct.dist)
## [1] 0.09650666

11.0.2 RH accordance rate

RHaccordance <- ddply(resultsRH.perTrial.all, c("condition", "subject"), summarise,
                          RUplus = length(which(resultsRH.perTrial.all$trialType == "1.5"))
                        )


RHaccordance <- resultsRH.perTrial.all %>% 
    group_by(subject, condition) %>%
    summarise(
      RUplus  = sum(trialType=="1.5"),
      RUminus = sum(trialType=="0.5"),
      RUrate = sum(trialType=="1.5")/(sum(trialType=="1.5") + sum(trialType=="0.5"))
      )
RHaccordance <- as.data.frame(RHaccordance)

mean(RHaccordance$RUrate[RHaccordance$condition == 1])
## [1] 0.8602743
mean(RHaccordance$RUrate[RHaccordance$condition == 2])
## [1] 0.6806216
sd(RHaccordance$RUrate[RHaccordance$condition == 1])
## [1] 0.1812772
sd(RHaccordance$RUrate[RHaccordance$condition == 2])
## [1] 0.2182993

11.0.3 α

mean(modelResults$alpha.1)
## [1] 0.7796724
mean(modelResults$alpha.2)
## [1] 0.4533456
sd(modelResults$alpha.1)
## [1] 0.0662246
sd(modelResults$alpha.2)
## [1] 0.1032093

11.0.4 RTs

# aggregate for each individual (1st level)
meanRTsindividual <- resultsRH.perTrial.all %>% 
    group_by(subject, condition, trialType) %>%
    summarise(
      RTs  = mean(RT)
      )
meanRTsindividual <- as.data.frame(meanRTsindividual)

# aggregate for each condition (2nd level)
meanRTs <- meanRTsindividual %>% 
    group_by(condition, trialType) %>%
    summarise(
      RTs  = mean(RTs)
      )
meanRTs <- as.data.frame(meanRTs)
meanRTs
##   condition trialType      RTs
## 1         1         0 2044.248
## 2         1       0.5 1882.800
## 3         1       1.5 1699.878
## 4         1         2 2236.263
## 5         2         0 2273.584
## 6         2       0.5 2412.643
## 7         2       1.5 2216.989
## 8         2         2 2524.360
stdRTs <- meanRTsindividual %>% 
    group_by(condition, trialType) %>%
    summarise(
      SDs  = sd(RTs)
      )
stdRTs <- as.data.frame(stdRTs)
stdRTs
##   condition trialType       SDs
## 1         1         0  730.6068
## 2         1       0.5  679.5547
## 3         1       1.5  547.7141
## 4         1         2  865.1021
## 5         2         0  940.7744
## 6         2       0.5 1123.3310
## 7         2       1.5  904.1902
## 8         2         2  938.0043

11.0.5 meanConfs

# 1st level
meanCONFsindividual <- resultsRH.perTrial.all %>% 
    group_by(subject, condition, trialType) %>%
    summarise(
      CONF  = mean(choiceConf)
      )
meanCONFsindividual <- as.data.frame(meanCONFsindividual)

# 2nd level
meanCONFs <- meanCONFsindividual %>% 
    group_by(condition, trialType) %>%
    summarise(
      meanCONFs  = mean(CONF)
      )
meanCONFs <- as.data.frame(meanCONFs)
meanCONFs
##   condition trialType meanCONFs
## 1         1         0  23.02080
## 2         1       0.5  32.05971
## 3         1       1.5  51.81117
## 4         1         2  54.23474
## 5         2         0  20.31187
## 6         2       0.5  28.38991
## 7         2       1.5  42.20778
## 8         2         2  52.12550
stdCONFs <- meanCONFsindividual %>% 
    group_by(condition, trialType) %>%
    summarise(
      stdCONFs  = sd(CONF)
      )
stdCONFs <- as.data.frame(stdCONFs)
stdCONFs
##   condition trialType stdCONFs
## 1         1         0 18.35739
## 2         1       0.5 22.78938
## 3         1       1.5 19.96125
## 4         1         2 17.57937
## 5         2         0 17.19012
## 6         2       0.5 18.92816
## 7         2       1.5 22.39717
## 8         2         2 22.88766
# Then run this from outside
# rmarkdown::render('/scripts/Filevichetal_2_SI_and_cleanSubjs.R')