TablesAndFigures.Rmd

---
title: "Tables and Figures"
output: html_document
---
Setup - install packages and source functions. 
```{r}
# clean slate to be sure
rm(list = ls())

# check for pacman package and install if not found
if (!require("pacman")) install.packages("pacman")
pacman::p_load(afex, readr, dplyr, tidyr, ggplot2, cowplot, devtools, here, tidyverse, magrittr, reshape2, rjags, coda, lattice, ggmcmc, foreach, doParallel, ggpubr, ggpol, patchwork, ggcorrplot, ggstatsplot, GGally, BayesFactor, gridExtra, viridis, gghalves, gt, rstatix, scales)

# load the metaSDT package from github
if (!require("metaSDT")) devtools::install_github("craddm/metaSDT")
# load Rousselet correlation package 
if (!require("bootcorci")) devtools::install_github("GRousselet/bootcorci")

source(here("r", "remove_outliers_robust.R"))
source(here("r", "outliers.R"))
source(here("r", "trials2counts.R"))
source(here("r", "metad_indiv.R"))
source(here("r", "fit_mle.R"))
source(here("r", "sdt_functions.R"))
source(here("r", "Function_metad_groupcorr.R"))
source(here("r", "BootCorci.R"))
source(here("r", "Function_scatterplot.R"))
source(here("r", "function_rhoPlot.R"))
source(here("r", "function_apa.R"))
source(here("r", "Bayesfactor.R"))

# set apa theme for plots
apatheme=theme_bw()+ #theme
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),
        panel.border=element_blank(),
        axis.line=element_line(),
        text = element_text(size = 15),
        axis.title = element_text(size = 12))



# spider-man color palette

colors = c("#DF1F2D", "#B11313", "#2B3784", "#447BBE")

```


Rearrange accuracy data for plotting summary variables

```{r}
 
#make the four bins: 
metadataPlot <- metadata %>%group_by(subject) %>% 
mutate(trial_number = row_number()) %>%
group_by(subject) #create a trial number column
  
#split in the four bins according to task 
metadataPlot <- metadataPlot %>%
    group_by(modality) %>%
    mutate(trial_bin = ntile(trial_number, n = 4)) 
  
#Summarise means for each bin and modality
metadataPlotMean <- metadataPlot %>%
group_by(trial_bin, modality) %>%
summarise(mean.acc = mean(accuracy),
          sd.acc = sd(accuracy), n.acc = n_distinct(subject)) %>%
mutate(se.acc = sd.acc / sqrt(n.acc),
         lower.ci.acc = mean.acc - qt(1 - (0.05 / 2), n.acc - 1) * se.acc,
         upper.ci.acc = mean.acc + qt(1 - (0.05 / 2), n.acc - 1) * se.acc) %>%
ungroup()


#Average confidence across domains, split by correct and incorrect trials
#split by correct and incorrect:
tmp_data <- metadata  %>%
   mutate(accuracy = factor(accuracy)) %>%
  mutate(accuracy = recode_factor(accuracy, '1' = "Hits", '0' = "Misses"))

#Average confidence for Correct trials:
ConfByError <- tmp_data  %>%
  group_by(subject, modality, accuracy) %>% 
  summarise(avgCF_Corr=mean(confidence, na_rm=TRUE)) 

```


Here we plot the Figures from the Paper: 
Figure 2: d' and group staircase: 
```{r}
# Figure2 

#Plot of d' without outliers: 

p1<-ggplot(metadata_exc, aes(x = modality, y = da, fill= modality)) +
   geom_half_point(range_scale = 8/10, shape = 21, alpha = 1/10, side = "l")+
   geom_half_boxplot(outlier.shape = NA, notch = TRUE, width = 1/4,
                     alpha = 1, color = "black", side = "r")+
  guides(fill="none") +
  ggtitle("Cognitive Sensitivity") +  labs(y = "D-prime", x = "Modality") + 
  theme_cowplot() + 
  scale_fill_manual(values = colors)+
  scale_x_discrete(labels = c("Calories", "GDP","Memory", "Vision")) +
  theme(legend.position = "none", 
        # aspect.ratio = 1,
        plot.title = element_text(size = 12), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=10,face="bold"))

#Plot of group-level Staircases: 


pd <- position_dodge(0.1)
p2 <- ggplot(metadataPlotMean, aes(x = trial_bin, y = mean.acc, group = modality, color= modality)) +
  geom_line(position = pd) +
  geom_point(position = pd) + 
  geom_errorbar(aes(ymin=lower.ci.acc,  ymax = upper.ci.acc),
                size = 1, width=0, position = pd, alpha = 1/4) +
  ylab("Average Accuracy")+
  xlab("Task Quarter") + 
  theme_cowplot() +
  ggtitle("Group-level Staircases") + 
  scale_linetype_discrete(name = "Modality", breaks = c("Calories", "GDP","Memory", "Vision")) +
  scale_x_continuous(labels = label_ordinal())+
  scale_y_continuous(breaks = seq(.6, .9, .1),
                     labels = scales::percent,
                     limits = c(.6, .9)) +
  theme(legend.title = element_blank(), 
       # aspect.ratio = 1,
        legend.position = c(1, 1), 
        legend.justification = c(1,1),
        legend.text = element_text(size = 10), 
        legend.direction = "horizontal", 
        plot.title = element_text(size = 12), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=10,face="bold")) +
  guides(colour = guide_legend(nrow = 2))+
  scale_colour_manual(values = colors)

f2 <- p1 + p2 

f2 <- f2 +  plot_annotation(tag_levels = 'A',   tag_suffix = ')') & 
  theme(plot.tag = element_text(size = 12))

f2
```


Figure 4
```{r}


p3 <- ggplot(data=ConfByError, 
       aes(x = modality, y = avgCF_Corr, fill=modality, shape = accuracy)) + 

  geom_half_point(aes(color = modality), range_scale = 1/2, alpha = 1/10,
                  side = "l", position = position_dodge(c(3/4,3/4)))+
  geom_half_boxplot(outlier.shape = NA, 
                    side = "r", notch = TRUE, width = 1/2,
                    alpha = 1, color = "black", show.legend = FALSE, 
                    position = position_dodge(c(3/4,3/4)))+
  theme_cowplot() +
  labs(y = "Average Confidence",
        x = "Modality", 
        fill = "Legend") + 
  scale_x_discrete(labels = c("Calories", "GDP","Memory", "Vision"))+ 
  scale_shape_manual(values=c(1, 4), guide = "legend")+
  guides(shape = guide_legend(override.aes = list(color = "black", alpha = 1))) +
  scale_fill_manual(values = colors, guide = "none")+
  scale_color_manual(values = colors, guide = "none")+
  scale_y_continuous(breaks = seq(1, 7, 1),
                     limits = c(.5, 7)) +
  ggtitle("Confidence by Accuracy") +
    theme(legend.position = c(.01,.05),
          legend.spacing.y = unit(1, "mm"),
          legend.spacing.x = unit(2, "mm"),
          legend.text.align = 0,
          legend.box.just = "center",
          #aspect.ratio = 1,
          legend.direction="horizontal",
          legend.background = element_blank(),
          legend.text = element_text(size = 10 ), 
          legend.title = element_blank(), 
          plot.title = element_text(size = 12), 
          axis.text=element_text(size=10),
          axis.title=element_text(size=10,face="bold")) 

p3
# M-ratio across modalities:


p4<-ggplot(metadata_exc, aes(x = modality, y = mratio, fill= modality)) +
    geom_hline(yintercept = 1, linetype="dashed", color = "black", size=.5)+
    aes(ymin = 0) +
    geom_blank() +
   geom_half_point(range_scale = 8/10, shape = 21, alpha = 1/10, side = "l")+
   geom_half_boxplot(outlier.shape = NA, notch = TRUE, width = 1/4, alpha = 1,
                     color = "black", side = "r")+
   guides(color="none") +
  ggtitle("Metacognitive Efficiency")+
  labs(y = "Meta-d'/d'", x = "Modality", color = "Legend") + 
  scale_fill_manual(values = colors)+
  scale_x_discrete(labels = c("Calories", "GDP","Memory", "Vision")) +
    scale_y_continuous(breaks = seq(-2, 4, 1),
                     limits = c(-2, 4)) +
   theme_cowplot() +
   theme(legend.position = "none", 
         #aspect.ratio = 1, 
         plot.title = element_text(size = 12), 
        axis.text=element_text(size=10),
        axis.title=element_text(size=10,face="bold")) 
  
  
 
  
f3 <- p3 + p4  +  plot_annotation(tag_levels = 'A',   tag_suffix = ')') & 
  theme(plot.tag = element_text(size = 10))

  f3
  


f3a <- (p1 + p2) / (p3 + p4) +  plot_annotation(tag_levels = 'A',   tag_suffix = ')') & 
  theme(plot.tag = element_text(size = 10))

f3a

ggsave(here("figs", "Figure 2.png"), height = 16.8, width = 16.8, units = "cm", dpi = 300)  

  
```




Figure 3:  
```{r}

#The scatter plot for the MLE estimation of Mratio: 
#Metacognitive Efficiency
#Indicate the estimated correlations: 

S1title <- paste("r =",round(res1.1$estimate, digits=3))
S2title <- paste("r =",round(res1.2$estimate, digits=3))
S3title <- paste("r =",round(res1.3$estimate, digits=3))
S4title <- paste("r =",round(res1.4$estimate, digits=3))
S5title <- paste("r =",round(res1.5$estimate, digits=3))
S6title <- paste("r =",round(res1.6$estimate, digits=3))


scat2 <- as_tibble(corr1)
s11 <- scatter(scat2, scat2$mratio_memory, scat2$mratio_Calories) + theme(plot.title = element_text( size = 10)) + scale_y_continuous("Calories")+scale_x_continuous("Memory") + ggtitle(paste(S1title))
s12 <- scatter(scat2, scat2$mratio_memory, scat2$mratio_vision) + theme(plot.title = element_text( size = 10)) + scale_y_continuous("Vision")+scale_x_continuous("Memory") + ggtitle(paste(S2title))
s13 <- scatter(scat2, scat2$mratio_memory, scat2$mratio_GDP) + theme(plot.title = element_text( size = 10)) + scale_y_continuous("GDP")+scale_x_continuous("Memory") + ggtitle(paste(S3title))
s14 <- scatter(scat2, scat2$mratio_GDP, scat2$mratio_vision) +  theme(plot.title = element_text( size = 10)) +scale_y_continuous("Vision")+scale_x_continuous("GDP") + ggtitle(paste(S4title))
s15 <- scatter(scat2, scat2$mratio_GDP, scat2$mratio_Calories) + theme(plot.title = element_text( size = 10)) + scale_y_continuous("Calories")+ scale_x_continuous("GDP") + ggtitle(paste(S5title))
s16 <- scatter(scat2, scat2$mratio_Calories, scat2$mratio_vision) + theme(plot.title = element_text( size = 10)) + scale_y_continuous("Vision")+scale_x_continuous("Calories") + ggtitle(paste(S6title))

comPlot2 <- ggarrange(s11, s12, s13, s14, s15, s16) 
annotate_figure(comPlot2, top = text_grob("Metacognitive Efficiency", color = "black", 
                                         face = "bold", size = 14))    

#The plots for the hierarhical model - sample estimations of rho: 

#For Astrid: 
output <- readRDS("C:/Users/Astrid/Documents/DomainGen/modelfit.RDS")

fit_filename = "/home/micah/metad_groupfit.RDS" ## change this to the correct place! 

if (file.exists(fit_filename) == 1 ) {

output <- readRDS(fit_filename)
  
} else{
  
message("OBS: You need to fit the hierarhical model using 'fit_hierarhicical_metamodel.rmd' - this will take a LONG time")

}

#Prepare the hierarhical fit for the figure: 
Value <- gelman.diag(output, confidence = 0.95)
Rhat <- data.frame(conv = Value$psrf)


# Values (mean and CI)
Value <- summary(output)
stat <- data.frame(mean = Value$statistics[,"Mean"])
stat %<>%
  rownames_to_column(var = "name") %>% 
  cbind(CILow = Value$quantiles[,"2.5%"]) %>% 
  cbind(CIUp = Value$quantiles[,"97.5%"])


# HDI function 
HDI <- data.frame(HPDinterval(output, prob = 0.95))
HDI %<>%
  rownames_to_column(var = "name")

# Posterior distributions 
mcmc.sample <- ggs(output) 


mcmc.rho <- mcmc.sample %>% 
  filter(Parameter == "rho[1]"| Parameter == "rho[2]"| Parameter == "rho[3]"| Parameter == "rho[4]"| Parameter == "rho[5]"| Parameter == "rho[6]")

#Plot titles: 
meanRho <- mcmc.rho %>% group_by(Parameter) %>% summarise(value= mean(value, na.rm=TRUE)) 

#Fit stats summary. 
Fit <- stat %>%
  cbind(lower = HDI$lower,
        upper = HDI$upper,
        Rhat = Rhat[,1])

H1title <- paste("Rho =",round(meanRho$value[1], digits=3))
H2title <- paste("Rho =",round(meanRho$value[2], digits=3))
H3title <- paste("Rho =",round(meanRho$value[3], digits=3))
H4title <- paste("Rho =",round(meanRho$value[4], digits=3))
H5title <- paste("Rho =",round(meanRho$value[5], digits=3))
H6title <- paste("Rho =",round(meanRho$value[6], digits=3))

# split per correlation pair: 
Rho1 <- mcmc.sample %>% filter(Parameter == "rho[1]") 
Rho2 <- mcmc.sample %>% filter(Parameter == "rho[2]") 
Rho3 <- mcmc.sample %>% filter(Parameter == "rho[3]") 
Rho4 <- mcmc.sample %>% filter(Parameter == "rho[4]") 
Rho5 <- mcmc.sample %>% filter(Parameter == "rho[5]") 
Rho6 <- mcmc.sample %>% filter(Parameter == "rho[6]") 

#plot Rho samples + indication of mean rho estimation 
R1 <- RhoPlot(Rho1, "rho[1]") + labs(title = "Memory and Calories", subtitle= paste(H1title))
R2 <- RhoPlot(Rho2, "rho[2]") + labs(title = "Memory and GDP", subtitle= paste(H2title))
R3 <- RhoPlot(Rho3, "rho[3]") + labs(title = "Memory and Vision", subtitle= paste(H3title))
R4 <- RhoPlot(Rho4, "rho[4]") + labs(title = "GDP and Calories", subtitle= paste(H4title))
R5 <- RhoPlot(Rho5, "rho[5]") + labs(title = "Calories and Vision", subtitle= paste(H5title))
R6 <- RhoPlot(Rho6, "rho[6]") + labs(title = "GDP and Vision", subtitle= paste(H6title))

#Combine plots for the hierarhical fit: 
RhoPlot <- ggarrange(R1, R3, R2, R6, R4, R5)
CorHi <- annotate_figure(RhoPlot, top = text_grob("Hierarhical Bayesian estimation", color = "black", face = "bold", size = 14)) 

#MLE fit:
CorMLE <- annotate_figure(comPlot2, top = text_grob("Maximum likelihood estimation  ", color = "black", face = "bold", size = 14)) 

#Combine MLe and hierarhical estimation of correlations: 
Figure3 <- ggarrange(CorMLE, CorHi, ncol= 1, nrow=2)  
  annotate_figure(Figure3, top = text_grob("Figure 3", color = "black",face = "bold", size = 14))

  
ggsave(here("figs", "Figure 3.png"), height = 7, width = 12)  

```

Figure 3 without Mratio outliers: 
```{r}

#exclude those where mratio could not be estimated (MLE):  
metadata_excMratio <- metadata_exc%>%filter(mratio!="NaN")

#mratio
metadata_excMratio <- metadata_excMratio%>%
  group_by(modality) %>%
  mutate(mratio = remove_outliers_robust(mratio)) 

#Pivot: 
fit_data_wideExc <- metadata_excMratio %>% pivot_wider(
  names_from = modality,
  values_from = c(da, mda, mratio, c, avg_conf)
)
  
#estimate correlation: 
corr1E <- fit_data_wideExc %>%
  na.omit() %>%
  select(starts_with("mratio"))
res1.1E <- corci(corr1E$mratio_memory, corr1E$mratio_Calories, method="spearman")
res1.2E <- corci(corr1E$mratio_memory, corr1E$mratio_vision, method="spearman")
res1.3E <- corci(corr1E$mratio_memory, corr1E$mratio_GDP, method="spearman")
res1.4E <- corci(corr1E$mratio_GDP, corr1E$mratio_vision, method="spearman")
res1.5E <- corci(corr1E$mratio_GDP, corr1E$mratio_Calories, method="spearman")
res1.6E <- corci(corr1E$mratio_Calories, corr1E$mratio_vision, method="spearman")

#put in dataframe: 
corrtableE <- data.frame(Metric=character(),
                 MemCalCoef=double(),
                 CI_lowerMC=double(),
                 CI_upperMC=double(),
                 P_valueMC=double(),
                 MemVisCoef=double(),
                 CI_lowerMV=double(),
                 CI_upperMV=double(),
                 P_valueMV=double(),
                 MemGDPCoef=double(),
                 CI_lowerMG=double(),
                 CI_upperMG=double(),
                 P_valueMG=double(),
                 GDPVisCoef=double(),
                 CI_lowerGS=double(),
                 CI_upperGS=double(),
                 P_valueGS=double(),
                 GDPCalCoef=double(),
                 CI_lowerGC=double(),
                 CI_upperGC=double(),
                 P_valueGC=double(),
                 CalVisCoef=double(),
                 CI_lowerCV=double(),
                 CI_upperCV=double(),
                 P_valueCV=double(),
                stringsAsFactors=FALSE)

corrtableE[1,] <- list("Mratio",res1.1E$estimate, res1.1E$conf.int[1], res1.1E$conf.int[2],
                      res1.1E$p.value, res1.2E$estimate,res1.2E$conf.int[1],
                      res1.2E$conf.int[2], res1.2E$p.value, res1.3E$estimate,
                      res1.3E$conf.int[1], res1.3E$conf.int[2], res1.3E$p.value,
                      res1.4E$estimate,res1.4E$conf.int[1], res1.4E$conf.int[2],
                      res1.4E$p.value, res1.5E$estimate,res1.5E$conf.int[1], 
                      res1.5E$conf.int[2], res1.5E$p.value, res1.6E$estimate,
                      res1.6E$conf.int[1],  res1.6E$conf.int[2], res1.6E$p.value)

S1Etitle <- paste("r =",round(res1.1E$estimate, digits=3))
S2Etitle <- paste("r =",round(res1.2E$estimate, digits=3))
S3Etitle <- paste("r =",round(res1.3E$estimate, digits=3))
S4Etitle <- paste("r =",round(res1.4E$estimate, digits=3))
S5Etitle <- paste("r =",round(res1.5E$estimate, digits=3))
S6Etitle <- paste("r =",round(res1.6E$estimate, digits=3))

#Plot the correlations: 
scat2E <- as_tibble(corr1E)
s11E <-scatter(scat2E, scat2E$mratio_memory, scat2E$mratio_Calories) + theme(plot.title = element_text( size = 10)) +
  scale_y_continuous("Calories")+scale_x_continuous("Memory") + ggtitle(paste(S1Etitle))
s12E <- scatter(scat2E, scat2E$mratio_memory, scat2E$mratio_vision) + theme(plot.title = element_text( size = 10)) +
  scale_y_continuous("Vision")+scale_x_continuous("Memory") + ggtitle(paste(S2Etitle))
s13E <- scatter(scat2E, scat2E$mratio_memory, scat2E$mratio_GDP) + theme(plot.title = element_text( size = 10)) +
  scale_y_continuous("GDP")+scale_x_continuous("Memory") + ggtitle(paste(S3Etitle))
s14E <- scatter(scat2E, scat2E$mratio_GDP, scat2E$mratio_vision) + theme(plot.title = element_text( size = 10)) +
  scale_y_continuous("Vision")+scale_x_continuous("GDP") + ggtitle(paste(S4Etitle))
s15E <- scatter(scat2E, scat2E$mratio_GDP, scat2E$mratio_Calories) + theme(plot.title = element_text( size = 10)) +
  scale_y_continuous("Calories")+ scale_x_continuous("GDP") + ggtitle(paste(S5Etitle))
s16E <- scatter(scat2E, scat2E$mratio_Calories, scat2E$mratio_vision) + theme(plot.title = element_text( size = 10)) +
  scale_y_continuous("Vision")+scale_x_continuous("Calories") + ggtitle(paste(S6Etitle))

comPlotE <- ggarrange(s11E, s12E, s13E, s14E, s15E, s16E) 
CorExcMLE<-annotate_figure(comPlotE, top = text_grob("Maximum likelihood estimation", color = "black", face = "bold", size = 14))   

#figure 3 now without MLE mratio outliers: 
Figure3E <- ggarrange(CorExcMLE, CorHi, ncol= 1, nrow=2)  
  annotate_figure(Figure3E, top = text_grob("Figure 3", color = "black",face = "bold", size = 14))

ggsave(here("figs", "Figure 3 without Outliers.png"), height = 7, width = 12)

```

Now the tables from the paper: 

Table 1: Demographics 
```{r}
#load the demographics information: 
survey_data_pre <- read_delim(here("data_summary", "/self_belief_pre_labels.csv"), 
                              ";", escape_double = FALSE, trim_ws = TRUE)

#Clean Data: 
survey_data_pre <- survey_data_pre %>% select(Sid, age, gender)

#Vector with IDs that should be included
ID_incl <- unique(metadata_exc$subject)

#demo only with included participants:
demo <- subset(survey_data_pre, is.element(survey_data_pre$Sid, ID_incl))

#summarise and format Gender:
demoG <- demo %>% group_by(gender) %>% 
  summarise(mean1=length(gender)) %>% mutate(per = 100*mean1/321) %>%
  mutate(Count=paste0(round(mean1), paste0(" (", round(per,2), " %)"))) %>% mutate(Gender=gender)%>%
  select(Gender, Count) 


#summarise and format age:
demoA <- demo %>% group_by(gender) %>% summarise(age1 =mean(age, na_rm=TRUE),
                         age_sd=sd(age, na.rm=TRUE)) %>% 
  mutate(Age= paste0(round(age1), paste0(" (", round(age_sd,3), " SD)"))) %>% select(Age, gender)


#MANUALLY CHANGE NON_BINARY AGE
demoA[3,1] <-'21 (0 SD)' 

#Join the two tables: 
table1 <- full_join(demoG, demoA, by = c("Gender"="gender"))
table1[3,1] <-'Non-binary'

#Apa format:
table1apa<-apa(table1, "Table 1: Demographics")
table1apa

```

Now the tables from the paper: 

Table 1: Demographics 
```{r}
#load the demographics information: 
survey_data_pre <- read_delim(here("data_summary", "/self_belief_pre_labels.csv"), 
                              ";", escape_double = FALSE, trim_ws = TRUE)

#Clean Data: 
survey_data_pre <- survey_data_pre %>% select(Sid, age, gender)

#Vector with IDs that should be included
ID_incl <- unique(metadata_exc$subject)

#demo only with included participants:
demo <- subset(survey_data_pre, is.element(survey_data_pre$Sid, ID_incl))

#summarise and format Gender:
demoG <- demo %>% group_by(gender) %>% 
  summarise(mean1=length(gender)) %>% mutate(per = 100*mean1/321) %>%
  mutate(Count=paste0(round(mean1), paste0(" (", round(per,2), " %)"))) %>% mutate(Gender=gender)%>%
  select(Gender, Count) 


#summarise and format age:
demoA <- demo %>% group_by(gender) %>% summarise(age1 =mean(age, na_rm=TRUE),
                         age_sd=sd(age, na.rm=TRUE)) %>% 
  mutate(Age= paste0(round(age1), paste0(" (", round(age_sd,3), " SD)"))) %>% select(Age, gender)


#MANUALLY CHANGE NON_BINARY AGE
demoA[3,1] <-'21 (0 SD)' 

#Join the two tables: 
table1 <- full_join(demoG, demoA, by = c("Gender"="gender"))
table1[3,1] <-'Non-binary'

#Apa format:
table1apa<-apa(table1, "Table 1: Demographics")
table1apa

```

Table 2: Hierarchical correlation of Metacognitive efficiency
NOTE you need the hierarhical fit to do this. 
```{r}

# select the HDIs for each rho  
rHDI <- HDI %>% filter(str_detect(name,"rho")) %>% select(name, lower, upper)
#join these HDIs with mean rho
RhoHDI <- full_join(meanRho, rHDI, 
                    by= c("Parameter"="name"))

#Change names and layout: 
Table2H <- RhoHDI %>% mutate(Pair= c("Memory & Calories", "Memory & Vision", "Memory & GDP", 
                                     "GDP & Vision", "GDP & Calories", "Calories & Vision"),
                             Rho=round(value,3),
                             HDI = paste(round(lower,3), round(upper,3), sep=";")
                             )

Table2H <- Table2H %>% mutate(`Rho (HDI)`= paste0(Rho, " (", HDI, ")" )) %>% select(Pair, `Rho (HDI)`)


#Put into APA format. 

Table2Hapa<-apa(Table2H, "Table 3: Hierarhical Bayesian estimation of correlation of Metacognitive Efficiency")

######Save tables. 
#######################
###

Table2Hapa
  
```

 MLE and Hierarchical fit of Mratio correlations:
```{r}
table23 
Table2H

TableHM <- full_join(table23, Table2H, by=c("Pair"="Pair"))
TableHM <- TableHM %>% mutate( `MLE fit` = paste0("r(314) = ", r, " ",`p-value (CI)`),
                               `Hierarchical fit` = `Rho (HDI)`) %>% 
  select(Pair, `MLE fit`, `Hierarchical fit`)

TableHMapa<-apa(TableHM, "Table 3: MLE and Hierarhical Bayesian estimation of correlation of Metacognitive Efficiency")

######Save tables. 
#######################
###

TableHMapa
```



Table 2: MLE correlation between domains 
```{r}

#combine upper and lower limit of CI into one. 
table2 <- corrtable %>% mutate(MC_CI = paste(round(corrtable$CI_lowerMC,3), round(corrtable$CI_upperMC,3), sep=";")) %>%
  mutate(MV_CI= paste(round(corrtable$CI_lowerMV,3), round(corrtable$CI_upperMV,3), sep=";")) %>%
  mutate(MG_CI= paste(round(corrtable$CI_lowerMG,3), round(corrtable$CI_upperMG,3), sep=";")) %>%
  mutate(GS_CI= paste(round(corrtable$CI_lowerGS,3), round(corrtable$CI_upperGS,3), sep=";")) %>%
  mutate(GC_CI= paste(round(corrtable$CI_lowerGC,3), round(corrtable$CI_upperGC,3), sep=";")) %>%
  mutate(CV_CI= paste(round(corrtable$CI_lowerCV,3), round(corrtable$CI_upperCV,3), sep=";")) 

#Select relevant columns 
table2 <- table2 %>% select(Metric, 
                            MemCalCoef, P_valueMC, MC_CI, 
                            MemVisCoef, P_valueMV, MV_CI, 
                            MemGDPCoef, P_valueMG, MG_CI, 
                            GDPVisCoef, P_valueGS, GS_CI,
                            GDPCalCoef, P_valueGC, GC_CI, 
                            CalVisCoef, P_valueCV, CV_CI
                            )

#Change the layout of the table: 

##Sensitivity 
#Table with correlation coefficents for all pairwise correlations: 
table22 <- data.frame(Pair=character(),
                 r=double(),
                 p=double(),
                 'Conf Interval'=double(),
                stringsAsFactors=FALSE)

#Insert vaules in dataframe: 
table22[1,] <- list("Memory & Calories",table2$MemCalCoef[1],table2$P_valueMC[1], table2$MC_CI[1])
table22[2,] <- list("Memory & Vision",table2$MemVisCoef[1],table2$P_valueMV[1], table2$MV_CI[1])
table22[3,] <- list("Memory & GDP",table2$MemGDPCoef[1],table2$P_valueMG[1], table2$MG_CI[1])
table22[4,] <- list("GDP & Vision",table2$GDPVisCoef[1],table2$P_valueGS[1], table2$GS_CI[1])
table22[5,] <- list("GDP & Calories",table2$GDPCalCoef[1],table2$P_valueGC[1], table2$GC_CI[1])
table22[6,] <- list("Calories & Vision",table2$CalVisCoef[1],table2$P_valueCV[1], table2$CV_CI[1])

#Combine p and C-value: 
table22 <- table22 %>% mutate(`p-value (CI)`= paste0('p = ', table22$p, ' (', table22$Conf.Interval, ')')) %>% select(Pair,r, `p-value (CI)`)




##Efficiency
#Table with correlation coefficents for all pairwise correlations: 
table23 <- data.frame(Pair=character(),
                 r=double(),
                 p=double(),
                 'Conf Interval'=double(),
                stringsAsFactors=FALSE)

#Insert vaules in dataframe: 
table23[1,] <- list("Memory & Calories",table2$MemCalCoef[2],table2$P_valueMC[2], table2$MC_CI[2])
table23[2,] <- list("Memory & Vision",table2$MemVisCoef[2],table2$P_valueMV[2], table2$MV_CI[2])
table23[3,] <- list("Memory & GDP",table2$MemGDPCoef[2],table2$P_valueMG[2], table2$MG_CI[2])
table23[4,] <- list("GDP & Vision",table2$GDPVisCoef[2],table2$P_valueGS[2], table2$GS_CI[2])
table23[5,] <- list("GDP & Calories",table2$GDPCalCoef[2],table2$P_valueGC[2], table2$GC_CI[2])
table23[6,] <- list("Calories & Vision",table2$CalVisCoef[2],table2$P_valueCV[2], table2$CV_CI[2])

#Combine p and C-value: 
table23 <- table23 %>% mutate(`p-value (CI)`= paste0('p = ', table23$p, ' (', table23$Conf.Interval, ')')) %>% select(Pair,r, `p-value (CI)`)

##Bias
#Table with correlation coefficents for all pairwise correlations: 
table24 <- data.frame(Pair=character(),
                 r=double(),
                 p=double(),
                 'Conf Interval'=double(),
                stringsAsFactors=FALSE)

#Insert vaules in dataframe: 
table24[1,] <- list("Memory & Calories",table2$MemCalCoef[3],table2$P_valueMC[3], table2$MC_CI[3])
table24[2,] <- list("Memory & Vision",table2$MemVisCoef[3],table2$P_valueMV[3], table2$MV_CI[3])
table24[3,] <- list("Memory & GDP",table2$MemGDPCoef[3],table2$P_valueMG[3], table2$MG_CI[3])
table24[4,] <- list("GDP & Vision",table2$GDPVisCoef[3],table2$P_valueGS[3], table2$GS_CI[3])
table24[5,] <- list("GDP & Calories",table2$GDPCalCoef[3],table2$P_valueGC[3], table2$GC_CI[3])
table24[6,] <- list("Calories & Vision",table2$CalVisCoef[3],table2$P_valueCV[3], table2$CV_CI[3])

#NOTE I MANUALLY CHANGE THE P-values to <.001 here! Thus with different data this may not be appropriate. 
table24 <- table24 %>% mutate(p = as.character("<.001"))
#Combine p and C-value: 
table24 <- table24 %>% mutate(`p-value (CI)`= paste0('p = ', table24$p, ' (', table24$Conf.Interval, ')')) %>% select(Pair,r, `p-value (CI)`)


#d' correlation: 
table25 <- data.frame(Pair=character(),
                 r=double(),
                 p=double(),
                 'Conf Interval'=double(),
                stringsAsFactors=FALSE)

#Insert vaules in dataframe: 
table25[1,] <- list("Memory & Calories",table2$MemCalCoef[4],table2$P_valueMC[4], table2$MC_CI[4])
table25[2,] <- list("Memory & Vision",table2$MemVisCoef[4],table2$P_valueMV[4], table2$MV_CI[4])
table25[3,] <- list("Memory & GDP",table2$MemGDPCoef[4],table2$P_valueMG[4], table2$MG_CI[4])
table25[4,] <- list("GDP & Vision",table2$GDPVisCoef[4],table2$P_valueGS[4], table2$GS_CI[4])
table25[5,] <- list("GDP & Calories",table2$GDPCalCoef[4],table2$P_valueGC[4], table2$GC_CI[4])
table25[6,] <- list("Calories & Vision",table2$CalVisCoef[4],table2$P_valueCV[4], table2$CV_CI[4])

#Combine p and C-value: 
table25 <- table25 %>% mutate(`p-value (CI)`= paste0('p = ', table25$p, ' (', table25$Conf.Interval, ')')) %>% select(Pair,r, `p-value (CI)`)
table25 <- table25 %>% mutate(r = round(r,3))
table25 <- table25 %>% mutate( `Spearman's Correlation` = paste0("r(314) = ", r, " ",`p-value (CI)`)) %>% select(Pair, `Spearman's Correlation`)

#Sort decimals: 
table22 <- table22 %>% mutate(r = round(r,3))
table23 <- table23 %>% mutate(r = round(r,3))
table24 <- table24 %>% mutate(r = round(r,3))


#Makes tables in APA format: 
CorEff<-apa(table23, "Table 2.a: MLE correlation of Metacognitive Efficiency")
CorSen<-apa(table22, "Table 2.b: MLE correlation of Metacognitive Sensitivity")
CorBias<-apa(table24, "Table 2.c: MLE correlation of Metacognitive Bias")
CorD <-apa(table25, "Table 2: MLE correlation of Cognitive Sensitivity")

######Save tables. 
#######################
###

CorBias
CorSen
CorEff
CorD

```


Table 3: Average cognitive and Metacognitive performance across tasks 
```{r}

#Average everything and calculate standard deviations: 
metadata_avg <- metadata_exc %>% 
  group_by(modality) %>% 
  summarise(dprime=mean(da, na_rm=TRUE), 
            criterion=mean(c, na_rm=TRUE),
            metad=mean(mda, na_rm=TRUE),
            mratio1=mean(mratio, na_rm=TRUE),
            avg_conf1=mean(avg_conf, na_rm=TRUE),
            dprime_sd=sd(da),
            criterion_sd=sd(c),
            metad_sd=sd(mda),
            mratio_sd=sd(mratio),
            avg_conf_sd=sd(avg_conf)
            )
    
#Change format: 
tableMLE <- metadata_avg %>% mutate(Modality=c("Calories", "GDP", "Memory", "Vision"),
                                  'dprime (SD)'=paste0((round(dprime,3)), paste0(" (", round(dprime_sd,3), ")")),
                                  'Criterion (SD)'=paste0((round(criterion,3)), paste0(" (", round(criterion_sd,3), ")")),
                                  'Meta-d (SD)'=paste0((round(metad,3)), paste0(" (", round(metad_sd,3), ")")),
                                  'Mratio (SD)'=paste0((round(mratio1,3)), paste0(" (", round(mratio_sd,3), ")")),
                                  'Avg Conf (SD)'=paste0((round(avg_conf1,3)), paste0(" (", round(avg_conf_sd,3), ")"))
                                  ) 
tableMLE <- tableMLE %>% select(Modality, 'dprime (SD)', 'Criterion (SD)', 'Avg Conf (SD)', 'Meta-d (SD)', 'Mratio (SD)')   

# Mratio estimations from Hierarhical Model.  

#Select Mratio and take the expontial of the value:
meanMratio <- Fit %>% filter(str_detect(name,"mu_logMratio")) %>%
  mutate(mean=exp(mean), 
         lower=exp(lower),
         upper=exp(upper))
  
#format and select: 
meanMratio <- meanMratio %>% mutate(Modality=c("Memory", "Calories", "GDP", "Vision"),
         'HMratio (HDI)'=paste0(round(mean,3), paste0(" (", round(lower,3),";",round(upper,3), ")"))) 
         
meanMratio <- meanMratio %>% select(Modality, 'HMratio (HDI)')
  
#Join the two tables MLE and Hierarhical dataframes. 

table3 <- left_join(tableMLE, meanMratio, 
                     by=c("Modality"="Modality"))


#Put into APA format.                                     
table3apa<-apa(table3, "Table 4: Cognitive and Metacognitive performance")

######Save tables. 
#######################
###

table3apa

```

Here are some additional plots - you can un-comment these sections to see them.
Note that the individual staircase plots may take a LONG time to produce. 

Plots of key metacognitive measures for both filtered and unfiltered data: 
```{r}
#Source here to get the additional plots of the distribution of d', criterion,
#meta-d', bias and m-ratio across condition, both before and after group-level
#exclusion: 

#source(here("r", "MLE_plots.R"))

#Plots showing the comparison of filtered and unfiltered plots: 
#annotate_figure(filunfil1, top = text_grob("Cognitive Criterion filtered/unfiltered", color = "black", face = "bold", size = 14))
#annotate_figure(filunfil5, top = text_grob("Cognitive Sensitivity filtered/unfiltered", color = "black",face = "bold", size = 14))
#annotate_figure(filunfil2, top = text_grob("Metaognitive Bias filtered/unfiltered", color = "black", face = "bold", size = 14))
#annotate_figure(filunfil3, top = text_grob("Metacognitive Sensitivity filtered/unfiltered",   color = "black", face = "bold", size = 14))
#annotate_figure(filunfil4, top = text_grob("Metacognitive Efficiency filtered/unfiltered", color = "black",  face = "bold", size = 14))

#Plot the unfiltered data´: 
#annotate_figure(Dataunfil1, top = text_grob("MLE Estimation of Cognitive performance (unfiltered)", color = "black",face = "bold", size = 14))  
#annotate_figure(Dataunfil2, top = text_grob("MLE Estimation of Metacognitive performance (unfiltered)", color = "black", face = "bold", size = 14))  

#Plot the filtered data´: 
#annotate_figure(Datafil1, top = text_grob("MLE Estimation of Cognitive performance (filtered)", color = "black",face = "bold", size = 14))  
#annotate_figure(Datafi2, top = text_grob("MLE Estimation of Metacognitive performance (unfiltered)", color = "black", face = "bold", size = 14)) 

```

Additional scatter plots and linear regressions: 
```{r}
#For scatters plots of MLE estimation of metacognitive bias and sensitvity, 
#source here: 

#source(here("r", "Scatterplots.R"))

#Plot them: 
#annotate_figure(comPlot1, top = text_grob("Metacognitive Sensitivity", color = "black", face = "bold", size = 14))  
#annotate_figure(comPlot3, top = text_grob("Metacognitive Bias", color = "black", face = "bold", size = 14))   

```


Individual-level staircase plots: (Takes a long time)
```{r}
# note that this procedure can take a LONG time 
#uncomment to make individual plots for staircase procedure:

#source(here("r", "make_indivi_plots.R"))

```