.Rhistory

FA2 <- FA2 / sum(FA2)
cum_H2 <- c(0, cumsum(H2))
cum_FA2 <- c(0, cumsum(FA2))
i <- 1
k <- numeric(Nratings)
for (c in 1:Nratings) {
k[i] <- (cum_H2[c + 1] - cum_FA2[c])^2 - (cum_H2[c] - cum_FA2[c + 1])^2
i <- i + 1
}
auroc2 <- 0.5 + 0.25 * sum(k)
# AUROC is stored in 'auroc2'
print(auroc2)
}
df$ResponseCorrect
df = df %>% drop_na()
correct = df$ResponseCorrect
conf = df$Confidence
Nratings = 4
i <- Nratings + 1
H2 <- numeric(Nratings)
FA2 <- numeric(Nratings)
for (c in 1:Nratings) {
H2[i-1] <- sum(conf == c & correct) + 0.5
FA2[i-1] <- sum(conf == c & !correct) + 0.5
i <- i - 1
}
H2 <- H2 / sum(H2)
FA2 <- FA2 / sum(FA2)
cum_H2 <- c(0, cumsum(H2))
cum_FA2 <- c(0, cumsum(FA2))
i <- 1
k <- numeric(Nratings)
for (c in 1:Nratings) {
k[i] <- (cum_H2[c + 1] - cum_FA2[c])^2 - (cum_H2[c] - cum_FA2[c + 1])^2
i <- i + 1
}
auroc2 <- 0.5 + 0.25 * sum(k)
# AUROC is stored in 'auroc2'
print(auroc2)
i
H2
H2 <- numeric(Nratings)
conf == c & correct
conf == c & !correct
c = 1
conf == c & !correct
conf == c & correct
conf = df$Confidence_bin
Nratings = 4
i <- Nratings + 1
H2 <- numeric(Nratings)
FA2 <- numeric(Nratings)
for (c in 1:Nratings) {
H2[i-1] <- sum(conf == c & correct) + 0.5
FA2[i-1] <- sum(conf == c & !correct) + 0.5
i <- i - 1
}
H2
FA2
df$Confidence_bin
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
modality == "Extero"
modality = "Extero"
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
correct = df$ResponseCorrect
conf = df$Confidence_bin
Nratings = 4
i <- Nratings + 1
H2 <- numeric(Nratings)
FA2 <- numeric(Nratings)
for (c in 1:Nratings) {
H2[i-1] <- sum(conf == c & correct) + 0.5
FA2[i-1] <- sum(conf == c & !correct) + 0.5
i <- i - 1
}
H2
FA2
H2 <- H2 / sum(H2)
FA2 <- FA2 / sum(FA2)
cum_H2 <- c(0, cumsum(H2))
cum_FA2 <- c(0, cumsum(FA2))
i <- 1
k <- numeric(Nratings)
for (c in 1:Nratings) {
k[i] <- (cum_H2[c + 1] - cum_FA2[c])^2 - (cum_H2[c] - cum_FA2[c + 1])^2
i <- i + 1
}
auroc2 <- 0.5 + 0.25 * sum(k)
# AUROC is stored in 'auroc2'
print(auroc2)
df
# AUROC is stored in 'auroc2'
print(auroc2)
m1 = glm(ResponseCorrect ~ Confidence_bin, data = df %>% filter(Modality == modality), family = binomial(link = "logit"))
roc_con = pROC::roc(df %>% filter(Modality == modality) %>% .$ResponseCorrect,m1$fitted)
auc = pROC::auc(roc_con)
# AUROC is stored in 'auroc2'
print(auroc2)
auc
pacman::p_load(tidyverse,ggdist,psycho,caret,patchwork, gt, cowplot, grid,reticulate,cmdstanr,posterior,rstan,bayesplot,here,rmarkdown)
np <- import("numpy")
#Here we read the same file as in the python notebook:
psychophysics_df = read_csv('https://github.com/embodied-computation-group/CardioceptionPaper/raw/main/data/Del2_merged.txt')
df = psychophysics_df %>% filter(Subject == "sub_0019")
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
df = psychophysics_df %>% filter(Subject == "sub_0019")
df$ResponseCorrect <- ifelse(df$ResponseCorrect == "", NA, df$ResponseCorrect)
df$Decision <- ifelse(df$Decision == "", NA, df$Decision)
df <- df %>%
mutate(
Decision = as.character(df$Decision),
ConfidenceRT = as.numeric(df$ConfidenceRT),
DecisionRT = as.numeric(df$DecisionRT),
Confidence = as.numeric(df$Confidence),
Condition = as.character(df$Condition),
listenBPM = as.numeric(df$listenBPM),
responseBPM = as.numeric(df$responseBPM),
ResponseCorrect = as.numeric(as.factor(df$ResponseCorrect)) - 1,
EstimatedThreshold = as.numeric(df$EstimatedThreshold),
EstimatedSlope = as.numeric(df$EstimatedSlope),
)
# check for NA's in the critical columns of the data:
# specify different ways of coding missing values:
missing_pres <- c(NA, "na", "N/A", "NaN", NaN, "n/a")
trials_missing <- df %>%
select(Decision, Confidence, ConfidenceRT, Confidence, Condition, listenBPM, responseBPM, nTrials) %>%
filter_all(any_vars(. %in% missing_pres)) %>%
.$nTrials
if (length(trials_missing) != 0) {
print(paste("Number of NA's = ", length(trials_missing), " detected in trials : "))
print(as.character(trials_missing))
}
# remove the NA's
df1 <- df %>% filter(!nTrials %in% trials_missing)
modality = "Intero"
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
df = df1
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
get_data = function(df, modality,bins){
if(bins == T){
m1 = glm(ResponseCorrect ~ Confidence_bin, data = df %>% filter(Modality == modality), family = binomial(link = "logit"))
}else{
m1 = glm(ResponseCorrect ~ Confidence, data = df %>% filter(Modality == modality), family = binomial(link = "logit"))
}
roc_con = pROC::roc(df %>% filter(Modality == modality) %>% .$ResponseCorrect,m1$fitted)
auc = pROC::auc(roc_con)
return(pROC::ggroc(roc_con)$data %>% mutate(Modality = modality, AUC = auc[[1]]))
}
model = get_data(df,modality, bins)
bins = T
model = get_data(df,modality, bins)
model
flemse(df$ResponseCorrect, df$Confidence_bin, 4)
model
0.8738095-0.8347739
df = df1
modality = "Extero"
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
model
0.8738095-0.8347739
flemse(df$ResponseCorrect, df$Confidence_bin, 4)
0.7892157-0.7382075
unique(psychophysics_df$Subject)
pacman::p_load(tidyverse,ggdist,psycho,caret,patchwork, gt, cowplot, grid,reticulate,cmdstanr,posterior,rstan,bayesplot,here,rmarkdown)
np <- import("numpy")
pacman::p_load(tidyverse,ggdist,psycho,caret,patchwork, gt, cowplot, grid,reticulate,cmdstanr,posterior,rstan,bayesplot,here,rmarkdown)
np <- import("numpy")
#Here we read the same file as in the python notebook:
psychophysics_df = read_csv('https://github.com/embodied-computation-group/CardioceptionPaper/raw/main/data/Del2_merged.txt')
df = psychophysics_df %>% filter(Subject == "sub_0175")
df$ResponseCorrect <- ifelse(df$ResponseCorrect == "", NA, df$ResponseCorrect)
df$Decision <- ifelse(df$Decision == "", NA, df$Decision)
df <- df %>%
mutate(
Decision = as.character(df$Decision),
ConfidenceRT = as.numeric(df$ConfidenceRT),
DecisionRT = as.numeric(df$DecisionRT),
Confidence = as.numeric(df$Confidence),
Condition = as.character(df$Condition),
listenBPM = as.numeric(df$listenBPM),
responseBPM = as.numeric(df$responseBPM),
ResponseCorrect = as.numeric(as.factor(df$ResponseCorrect)) - 1,
EstimatedThreshold = as.numeric(df$EstimatedThreshold),
EstimatedSlope = as.numeric(df$EstimatedSlope),
)
# check for NA's in the critical columns of the data:
# specify different ways of coding missing values:
missing_pres <- c(NA, "na", "N/A", "NaN", NaN, "n/a")
trials_missing <- df %>%
select(Decision, Confidence, ConfidenceRT, Confidence, Condition, listenBPM, responseBPM, nTrials) %>%
filter_all(any_vars(. %in% missing_pres)) %>%
.$nTrials
if (length(trials_missing) != 0) {
print(paste("Number of NA's = ", length(trials_missing), " detected in trials : "))
print(as.character(trials_missing))
}
# remove the NA's
df1 <- df %>% filter(!nTrials %in% trials_missing)
modality = "Extero"
get_data = function(df, modality,bins){
if(bins == T){
m1 = glm(ResponseCorrect ~ Confidence_bin, data = df %>% filter(Modality == modality), family = binomial(link = "logit"))
}else{
m1 = glm(ResponseCorrect ~ Confidence, data = df %>% filter(Modality == modality), family = binomial(link = "logit"))
}
roc_con = pROC::roc(df %>% filter(Modality == modality) %>% .$ResponseCorrect,m1$fitted)
auc = pROC::auc(roc_con)
return(pROC::ggroc(roc_con)$data %>% mutate(Modality = modality, AUC = auc[[1]]))
}
df = df1
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
bins = T
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
model
flemse(correct = df$ResponseCorrect, conf = df$Confidence_bin, Nratings = 4)
0.8104575-0.7658662
modality
flemse(correct = df$ResponseCorrect, conf = df$Confidence, Nratings = unique(df$Confidence))
df = df1
modality = "Intero"
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
flemse(correct = df$ResponseCorrect, conf = df$Confidence, Nratings = 4)
model
correct = df$ResponseCorrect
flemse(correct = df$ResponseCorrect, conf = df$Confidence_bin, Nratings = 4)
conf = df$Confidence_bin
Nratings = 4
correct = df$ResponseCorrect
i <- Nratings + 1
H2 <- numeric(Nratings)
FA2 <- numeric(Nratings)
H2
flemse(correct = df$ResponseCorrect, conf = df$Confidence_bin, Nratings = 4)
print(model$AUC[1])
flemse(correct = df$ResponseCorrect, conf = df$Confidence_bin, Nratings = 4)
pacman::p_load(tidyverse,ggdist,psycho,caret,patchwork, gt, cowplot, grid,reticulate,cmdstanr,posterior,rstan,bayesplot,here,rmarkdown)
np <- import("numpy")
#Here we read the same file as in the python notebook:
psychophysics_df = read_csv('https://github.com/embodied-computation-group/CardioceptionPaper/raw/main/data/Del2_merged.txt')
df = psychophysics_df %>% filter(Subject == "sub_0202")
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
results = single_sub_analysis(df,                                                                  #The raw dataframe
interoPost = interoPost,                                             #numpy array for the intero (NA if not avaliable)
exteroPost = exteroPost,                                             #numpy array for the extero (NA if not avaliable)
bayesian = F,                                                        #Bayesian Analysis (TRUE/FALSE)
model = NA,                                                          #Bayesian model here a stan script (NA if Bayesian is FALSE)
out = here::here("docs","source","examples","R"))                    #Output directory for results
#This line reads in a subject result file:
df = read_csv("https://raw.githubusercontent.com/embodied-computation-group/Cardioception/master/docs/source/examples/templates/data/HRD/HRD_final.txt")
#These next lines read in the psi posterior distributions for each modality.
#These will be saved with names depending on your conditions, so modify as needed.
interoPost = np$load(here("docs","source","examples","templates","data","HRD","Intero_posterior.npy"))
exteroPost = np$load(here("docs","source","examples","templates","data","HRD","Extero_posterior.npy"))
results = single_sub_analysis(df,                                                                  #The raw dataframe
interoPost = interoPost,                                             #numpy array for the intero (NA if not avaliable)
exteroPost = exteroPost,                                             #numpy array for the extero (NA if not avaliable)
bayesian = F,                                                        #Bayesian Analysis (TRUE/FALSE)
model = NA,                                                          #Bayesian model here a stan script (NA if Bayesian is FALSE)
out = here::here("docs","source","examples","R"))                    #Output directory for results
pacman::p_load(tidyverse,ggdist,psycho,caret,patchwork, gt, cowplot, grid,reticulate, here,rmarkdown)
np <- import("numpy")
#This line reads in a subject result file:
df = read_csv("https://raw.githubusercontent.com/embodied-computation-group/Cardioception/master/docs/source/examples/templates/data/HRD/HRD_final.txt")
#These next lines read in the psi posterior distributions for each modality.
#These will be saved with names depending on your conditions, so modify as needed.
interoPost = np$load(here("docs","source","examples","templates","data","HRD","Intero_posterior.npy"))
exteroPost = np$load(here("docs","source","examples","templates","data","HRD","Extero_posterior.npy"))
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
results = single_sub_analysis(df,                                                                  #The raw dataframe
interoPost = interoPost,                                             #numpy array for the intero (NA if not avaliable)
exteroPost = exteroPost,                                             #numpy array for the extero (NA if not avaliable)
bayesian = F,                                                        #Bayesian Analysis (TRUE/FALSE)
model = NA,                                                          #Bayesian model here a stan script (NA if Bayesian is FALSE)
out = here::here("docs","source","examples","R"))                    #Output directory for results
pacman::p_load(tidyverse,ggdist,psycho,caret,patchwork, gt, cowplot, grid,reticulate,cmdstanr,posterior,rstan,bayesplot,here,rmarkdown)
np <- import("numpy")
#Here we read the same file as in the python notebook:
psychophysics_df = read_csv('https://github.com/embodied-computation-group/CardioceptionPaper/raw/main/data/Del2_merged.txt')
df = psychophysics_df %>% filter(Subject == "sub_0202")
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
unique(psychophysics_df$Subject)
df = psychophysics_df %>% filter(Subject == "sub_0172")
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
results
df = psychophysics_df %>% filter(Subject == "sub_0067")
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
model = cmdstan_model(here("docs","source","examples","R","src","first_model.stan"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
results$AUC_plot
df = psychophysics_df %>% filter(Subject == "sub_0098")
model = cmdstan_model(here("docs","source","examples","R","src","first_model.stan"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
df = psychophysics_df %>% filter(Subject == "sub_0044")
model = cmdstan_model(here("docs","source","examples","R","src","first_model.stan"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
df = psychophysics_df %>% filter(Subject == "sub_0042")
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
results$AUC_plot
results$histogram_plot
results$conf_plot
df$ResponseCorrect <- ifelse(df$ResponseCorrect == "", NA, df$ResponseCorrect)
df$Decision <- ifelse(df$Decision == "", NA, df$Decision)
df <- df %>%
mutate(
Decision = as.character(df$Decision),
ConfidenceRT = as.numeric(df$ConfidenceRT),
DecisionRT = as.numeric(df$DecisionRT),
Confidence = as.numeric(df$Confidence),
Condition = as.character(df$Condition),
listenBPM = as.numeric(df$listenBPM),
responseBPM = as.numeric(df$responseBPM),
ResponseCorrect = as.numeric(as.factor(df$ResponseCorrect)) - 1,
EstimatedThreshold = as.numeric(df$EstimatedThreshold),
EstimatedSlope = as.numeric(df$EstimatedSlope),
)
# check for NA's in the critical columns of the data:
# specify different ways of coding missing values:
missing_pres <- c(NA, "na", "N/A", "NaN", NaN, "n/a")
trials_missing <- df %>%
select(Decision, Confidence, ConfidenceRT, Confidence, Condition, listenBPM, responseBPM, nTrials) %>%
filter_all(any_vars(. %in% missing_pres)) %>%
.$nTrials
if (length(trials_missing) != 0) {
print(paste("Number of NA's = ", length(trials_missing), " detected in trials : "))
print(as.character(trials_missing))
}
# remove the NA's
df1 <- df %>% filter(!nTrials %in% trials_missing)
df = df1
modality = "Extero"
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
get_data = function(df, modality,bins){
if(bins == T){
m1 = glm(ResponseCorrect ~ Confidence_bin, data = df %>% filter(Modality == modality), family = binomial(link = "logit"))
}else{
m1 = glm(ResponseCorrect ~ Confidence, data = df %>% filter(Modality == modality), family = binomial(link = "logit"))
}
roc_con = pROC::roc(df %>% filter(Modality == modality) %>% .$ResponseCorrect , m1$fitted)
auc = pROC::auc(roc_con)
return(pROC::ggroc(roc_con)$data %>% mutate(Modality = modality, AUC = auc[[1]]))
}
bins = T
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
print(model$AUC[1])
flemse(correct = df$ResponseCorrect, conf = df$Confidence_bin, Nratings = 4)
df$ResponseCorrect
df$Confidence_bin
df
nrow(df)
df <- df %>%
filter(Modality == modality) %>%
mutate(stimuli = as.numeric(Alpha > 0), Modality = as.factor(Modality), reponse = as.factor(as.numeric((Decision == "More")))) %>%
ungroup()
df$Confidence_bin <- discretebins(df, 4)[[1]]
model = get_data(df,modality, bins)
print(model$AUC[1])
flemse(correct = df$ResponseCorrect, conf = df$Confidence_bin, Nratings = 4)
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
model = cmdstan_model(here("docs","source","examples","R","src","first_model.stan"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
pacman::p_load(tidyverse,ggdist,psycho,caret,patchwork, gt, cowplot, grid,reticulate,cmdstanr,posterior,rstan,bayesplot,here,rmarkdown)
np <- import("numpy")
#Here we read the same file as in the python notebook:
psychophysics_df = read_csv('https://github.com/embodied-computation-group/CardioceptionPaper/raw/main/data/Del2_merged.txt')
df = psychophysics_df %>% filter(Subject == "sub_0042")
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
model = cmdstan_model(here("docs","source","examples","R","src","first_model.stan"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
#loading the functions to do the analysis:
source(here("docs","source","examples","R","src","firstlevelanalysis.R"))
model = cmdstan_model(here("docs","source","examples","R","src","first_model.stan"))
results = single_sub_analysis(df,
interoPost = NA,
exteroPost = NA,
bayesian = F,
model = model,
out = here::here("docs","source","examples","R"))
results$conf_plot
results$AUC_plot