TdcsFrameBetween/Behavioral data + R-script/R-Skript_FrameBetween_publication.R

#### FrameBetween ####


# # Sind alle notwendigen Pakete installiert?
# install.packages("car")
# install.packages("ez")
# install.packages("ggplot2")
# install.packages("pastecs")
# install.packages("tidyr")
# install.packages("Hmisc")
# install.packages("dplyr")
# install.packages("tidyverse")
# install.packages("gdata")
# install.packages("lme4")
# install.packages("psych")
# install.packages("psychReport")
# install.packages("rcompanion")
# install.packages("data.table")
# install.packages("hrbrthemes")
# install.packages("viridis")
library(ggplot2)
library(psychReport)
library(dplyr)
library(apa)
library(afex)
library(dplyr)
library(forcats)
library(data.table)
library(psych)
library(doBy)
library(rcompanion)
library(effectsize)



#### Read data ####


setwd("Your Directory")    

list_of_files <- list.files(path = "Your Directory", recursive = TRUE,
                            pattern = "\\.dat$", 
                            full.names = TRUE)

FrameBetween <- vroom(list_of_files)

FrameBetween <- subset(FrameBetween, FrameBetween$subjID !=2) # MEG artifacts

# stim: 1 = anodal /excitatory    2 = Sham/Placebo
FrameBetween$stim <- ifelse(FrameBetween$stim == 1, "Excitatory", "Sham")


#### Check of Programming ####
table  (FrameBetween$InitRew)
table  (FrameBetween$RisktoLoose)
table  (FrameBetween$GainLossTr)
table  (FrameBetween$subjID)
table  (FrameBetween$choice)
describe(FrameBetween$choice)
table  (FrameBetween$stim)
kreuztabelle <- xtabs (~ FrameBetween$RisktoLoose + FrameBetween$GainLossTr)
xtabs (~ FrameBetween$stim + FrameBetween$subjID)
print(kreuztabelle)

table(FrameBetween$subjID)

# Half of the trials risky choices?
observed_proportion <- c(10808, 12232)
theoretical_proportion <- c(1/2, 1/2)
chisq.test(observed_proportion, p=theoretical_proportion)

FrameBetween = FrameBetween %>%
  dplyr::group_by(subjID) %>%
  dplyr::mutate(TrialNumber = (row_number()-1) %% 320 +1) %>%
  ungroup()




#### Add expected value ####
FrameBetween <- FrameBetween %>% 
  mutate(expectation_valueRisk = case_when(RisktoLoose == 20 & InitRew == 100 ~ 80,
                                           RisktoLoose == 40 & InitRew == 100 ~ 60,
                                           RisktoLoose == 60 & InitRew == 100 ~ 40,
                                           RisktoLoose == 80 & InitRew == 100 ~ 20,
                                           RisktoLoose == 20 & InitRew == 75 ~ 60,
                                           RisktoLoose == 40 & InitRew == 75 ~ 45,
                                           RisktoLoose == 60 & InitRew == 75 ~ 30,
                                           RisktoLoose == 80 & InitRew == 75 ~ 15,
                                           RisktoLoose == 20 & InitRew == 50 ~ 40,
                                           RisktoLoose == 40 & InitRew == 50 ~ 30,
                                           RisktoLoose == 60 & InitRew == 50 ~ 20,
                                           RisktoLoose == 80 & InitRew == 50 ~ 10,
                                           RisktoLoose == 20 & InitRew == 25 ~ 20,
                                           RisktoLoose == 40 & InitRew == 25 ~ 15,
                                           RisktoLoose == 60 & InitRew == 25 ~ 10,
                                           RisktoLoose == 80 & InitRew == 25 ~ 5,
  ))

FrameBetween <- FrameBetween %>% 
  mutate(expectation_valueSafe = case_when(InitRew == 100 ~ 40,
                                           InitRew == 75 ~ 30,
                                           InitRew == 50 ~ 20,
                                           InitRew == 25 ~ 10,
  ))

FrameBetween <- FrameBetween %>% 
  mutate(expectation_valueChoice = case_when(choice == 1 ~ expectation_valueRisk,
                                             choice == 0 ~ expectation_valueSafe,
  ))


FrameBetween <- FrameBetween %>% 
  mutate(RationalDecision = case_when(choice == 1 & expectation_valueRisk >= expectation_valueSafe  ~ 1,
                                      choice == 0 & expectation_valueSafe >= expectation_valueRisk ~ 1,
                                      choice == 1 & expectation_valueRisk < expectation_valueSafe  ~ 0,
                                      choice == 0 & expectation_valueSafe < expectation_valueRisk ~ 0,
  ))


FrameBetween <- FrameBetween %>% 
  mutate(Outcome = case_when(choice == 1 & RiskFeedback == 22 | choice == 1 &  RiskFeedback == 52 | choice == 1 & RiskFeedback == 72 | choice == 1 & RiskFeedback == 102  ~ 1,
                             choice == 1 & RiskFeedback == 21 | choice == 1 & RiskFeedback == 51 | choice == 1 & RiskFeedback == 71 | choice == 1 & RiskFeedback == 101  ~ 0,
                             choice == 0 & GainLossTr == 0 & RiskFeedback == 0  ~ 1, # 99 oder 1
                             choice == 0 & GainLossTr == 1 & RiskFeedback == 0 ~ 0 # 99 oder 0
  ))

FrameBetween$choicepercent <- FrameBetween$choice*100

#FrameBetween <- FrameBetween %>% 
#  mutate(Outcome4 = case_when(choice == 1 & RiskFeedback == 22 | choice == 1 &  RiskFeedback == 52 | choice == 1 & RiskFeedback == 72 | choice == 1 & RiskFeedback == 102  ~ 4,
#                             choice == 1 & RiskFeedback == 21 | choice == 1 & RiskFeedback == 51 | choice == 1 & RiskFeedback == 71 | choice == 1 & RiskFeedback == 101  ~ 1,
#                             choice == 0 & GainLossTr == 0 & RiskFeedback == 0  ~ 3, 
#                             choice == 0 & GainLossTr == 1 & RiskFeedback == 0 ~ 2 
#  ))

FrameBetween <- data.table(FrameBetween)
FrameBetween[ , previousTr := shift(Outcome)]
#FrameBetween[ , previousTr4 := shift(Outcome4)]

PrevTrial <- subset(FrameBetween, FrameBetween$previousTr !=99)
PrevTrial <- subset(PrevTrial, PrevTrial$previousTr !="NA")
PrevTrial <- subset(PrevTrial, PrevTrial$TrialNumber !=1)
#PrevTrial4 <- subset(FrameBetween, FrameBetween$previousTr4 !=99)
#PrevTrial4 <- subset(PrevTrial, PrevTrial$previousTr4 !="NA")
PrevTrial$previousTr <- ifelse(PrevTrial$previousTr == 1, "Gain", "Loss")


# Choice = 0 --> no risk, Choice = 1 --> risk

# GainLossTr: 0 = GainTrial    1 = LossTrial

FrameRat <- subset(FrameBetween, FrameBetween$RisktoLoose !=60)

GLMRat <- glm(formula = RationalDecision ~ stim, family = binomial(logit), data = FrameRat)
summary(GLMRat)

GLMFrame <- glm(formula = choice ~ GainLossTr, family = binomial(logit), data = FrameBetween)
summary(GLMFrame)
effectsizes <- standardize_parameters(GLMFrame, exp = TRUE)

GLMIntFrame <- glm(formula = choice ~ stim * GainLossTr, family = binomial(logit), data = FrameBetween)
summary(GLMIntFrame)
effectsizes <- standardize_parameters(GLMIntFrame, exp = TRUE)

GLMALL <- glm(formula = choice ~ stim * GainLossTr * RisktoLoose, family = binomial(logit), data = FrameBetween)
summary(GLMALL)
effectsizes <- standardize_parameters(GLMALL, exp = TRUE)


GLMInitRew <- glm(formula = choice ~ stim * InitRew, family = binomial(logit), data = FrameBetween)
summary(GLMInitRew)


Frame_Sham <- subset(FrameBetween, FrameBetween$stim !="Excitatory")
Frame_Ano <- subset(FrameBetween, FrameBetween$stim !="Sham")

GLMInitRew <- glm(formula = choice ~ InitRew, family = binomial(logit), data = Frame_Sham)
summary(GLMInitRew)
effectsizes <- standardize_parameters(GLMInitRew, exp = TRUE)

GLMInitRew <- glm(formula = choice ~ InitRew, family = binomial(logit), data = Frame_Ano)
summary(GLMInitRew)

GLMprevTr <- glm(formula = choice ~ stim * previousTr, family = binomial(logit), data = PrevTrial)
summary(GLMprevTr)
effectsizes <- standardize_parameters(GLMprevTr, exp = TRUE)

Frame_Sham <- subset(PrevTrial, PrevTrial$stim !="Excitatory")
Frame_Ano <- subset(PrevTrial, PrevTrial$stim !="Sham")

Frame_prevLoss <- subset(PrevTrial, PrevTrial$previousTr !="Gain")
Frame_prevGain <- subset(PrevTrial, PrevTrial$previousTr !="Loss")


GLMprevTr <- glm(formula = choice ~ previousTr, family = binomial(logit), data = Frame_Ano)
summary(GLMprevTr)
effectsizes <- standardize_parameters(GLMprevTr, exp = TRUE)

GLMprevTr <- glm(formula = choice ~ previousTr, family = binomial(logit), data = Frame_Sham)
summary(GLMprevTr)
effectsizes <- standardize_parameters(GLMprevTr, exp = TRUE)

GLMprevTr <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame_prevLoss)
summary(GLMprevTr)

GLMprevTr <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame_prevGain)
summary(GLMprevTr)

GLMNull <- glm(formula = choice ~ 1, family = binomial(logit), data = FrameBetween)
summary(GLMNull)

anova(GLMFrame, GLMNull, test = "Chisq")
anova(GLMRat, GLMNull, test = "Chisq")
anova(GLMALL, GLMNull, test = "Chisq")



#### Post-Hocs Interaction Risk and Expected Value ####


apa(t.test(expectation_valueChoice ~ stim, data = FrameBetween, var.equal = TRUE, paired = FALSE,
           alternative = "greater")) # two.sided

groupwiseMean(expectation_valueChoice ~ stim,
              data = FrameBetween,
              conf = 0.95,
              digits = 3)

FrameAgg <- aggregate(choice ~ subjID + stim + RisktoLoose,
                      data = FrameBetween, FUN = mean)


Frame20 <- subset(FrameBetween, FrameBetween$RisktoLoose !=80)
Frame20 <- subset(Frame20, Frame20$RisktoLoose !=60)
Frame20 <- subset(Frame20, Frame20$RisktoLoose !=40)

apa(t.test(expectation_valueChoice ~ stim, data = Frame20, var.equal = TRUE, paired = FALSE,
                       alternative = "greater")) # two.sided

GLM20 <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame20)
summary(GLM20)
effectsizes <- standardize_parameters(GLM20, exp = TRUE)

xtabs (~ Frame20$choice + Frame20$stim)

prop.test(x = c(1089, 1001), n = c(1280, 1280),alternative = "two.sided")

Frame40 <- subset(FrameBetween, FrameBetween$RisktoLoose !=80)
Frame40 <- subset(Frame40, Frame40$RisktoLoose !=60)
Frame40 <- subset(Frame40, Frame40$RisktoLoose !=20)

apa(t.test(expectation_valueChoice ~ stim, data = Frame40, var.equal = TRUE, paired = FALSE,
           alternative = "greater")) # two.sided

xtabs (~ Frame40$choice + Frame40$stim)

GLM40 <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame40)
summary(GLM40)

prop.test(x = c(852, 867), n = c(1280, 1280),alternative = "two.sided")

Frame60 <- subset(FrameBetween, FrameBetween$RisktoLoose !=80)
Frame60 <- subset(Frame60, Frame60$RisktoLoose !=20)
Frame60 <- subset(Frame60, Frame60$RisktoLoose !=40)

apa(t.test(expectation_valueChoice ~ stim, data = Frame60, var.equal = TRUE, paired = FALSE,
           alternative = "greater"))

xtabs (~ Frame60$choice + Frame60$stim)

GLM60 <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame60)
summary(GLM60)
effectsizes <- standardize_parameters(GLM60, exp = TRUE)

prop.test(x = c(468, 637), n = c(1280, 1280),alternative = "two.sided")

Frame80 <- subset(FrameBetween, FrameBetween$RisktoLoose !=60)
Frame80 <- subset(Frame80, Frame80$RisktoLoose !=20)
Frame80 <- subset(Frame80, Frame80$RisktoLoose !=40)

apa(t.test(expectation_valueChoice ~ stim, data = Frame80, var.equal = TRUE, paired = FALSE,
           alternative = "greater"))

xtabs (~ Frame80$choice + Frame80$stim)

GLM80 <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame80)
summary(GLM80)
effectsizes <- standardize_parameters(GLM80, exp = TRUE)
interpret_oddsratio(1.88, rules = "cohen1988")

prop.test(x = c(223, 361), n = c(1280, 1280),alternative = "two.sided")

apa(t.test(expectation_valueChoice ~ stim, data = FrameRat, var.equal = TRUE, paired = FALSE,
           alternative = "greater")) # two.sided

#### Post-hoc Tests Interaction Frame ####
LossFrame <- subset(FrameBetween, FrameBetween$GainLossTr !="GainTr")
GainFrame <- subset(FrameBetween, FrameBetween$GainLossTr !="LossTr")

GainFrame$Diff <- GainFrame$choice - LossFrame$choice

bar <- ggplot(GainFrame, aes(stim, Diff, fill = stim))
bar + 
  stat_summary(fun = mean, geom = "bar", position = "dodge") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position = position_dodge(width = 0.90),width = 0.2) +
  labs (x = "Frame", y = "Risky Choice", fill = "Stim") +
  scale_fill_brewer(palette = "Set1")
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
      panel.background = element_blank(), axis.line = element_line(colour = "black"))

xtabs (~ GainFrame$choice + GainFrame$stim)

prop.test(x = c(1293, 1384), n = c(2560, 2560),alternative = "two.sided")

xtabs (~ LossFrame$choice + LossFrame$stim)

prop.test(x = c(1339, 1482), n = c(2560, 2560),alternative = "two.sided")

GainFrameexc <- subset(GainFrame, GainFrame$stim !="Sham")
GainFrameSham <- subset(GainFrame, GainFrame$stim !="Excitatory")

table(GainFrameexc$Diff)
table(GainFrameSham$Diff)

prop.test(x = c(640, 680), n = c(2560, 2560),alternative = "less")

#### Post-Hoc Interaction InitRew ####

Frame25 <- subset(FrameBetween, FrameBetween$InitRew !=50)
Frame25 <- subset(Frame25, Frame25$InitRew !=75)
Frame25 <- subset(Frame25, Frame25$InitRew !=100)

xtabs (~ Frame25$choice + Frame25$stim)

prop.test(x = c(652, 647), n = c(1280, 1280),alternative = "two.sided")

Frame50 <- subset(FrameBetween, FrameBetween$InitRew !=25)
Frame50 <- subset(Frame50, Frame50$InitRew !=75)
Frame50 <- subset(Frame50, Frame50$InitRew !=100)

xtabs (~ Frame50$choice + Frame50$stim)

prop.test(x = c(658, 692), n = c(1280, 1280),alternative = "two.sided")

Frame75 <- subset(FrameBetween, FrameBetween$InitRew !=50)
Frame75 <- subset(Frame75, Frame75$InitRew !=25)
Frame75 <- subset(Frame75, Frame75$InitRew !=100)

xtabs (~ Frame75$choice + Frame75$stim)

prop.test(x = c(652, 728), n = c(1280, 1280),alternative = "two.sided")

GLM75 <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame75)
summary(GLM75)
effectsizes <- standardize_parameters(GLM75, exp = TRUE)

Frame100 <- subset(FrameBetween, FrameBetween$InitRew !=25)
Frame100 <- subset(Frame100, Frame100$InitRew !=50)
Frame100 <- subset(Frame100, Frame100$InitRew !=75)

xtabs (~ Frame100$choice + Frame100$stim)

prop.test(x = c(670, 799), n = c(1280, 1280),alternative = "two.sided")

GLM100 <- glm(formula = choice ~ stim, family = binomial(logit), data = Frame100)
summary(GLM100)
effectsizes <- standardize_parameters(GLM100, exp = TRUE)


#### Post-Hoc previous Trial ####
prvLoss <- subset(PrevTrial, PrevTrial$previousTr !="Gain")
prvGain <- subset(PrevTrial, PrevTrial$previousTr !="Loss")
Exc <- subset(PrevTrial, PrevTrial$stim !="Sham")
Sham <- subset(PrevTrial, PrevTrial$stim !="Excitatory")

xtabs (~ prvLoss$choice + prvLoss$stim)

prop.test(x = c(1126, 1262), n = c(2217, 2181),alternative = "two.sided")

xtabs (~ prvGain$choice + prvGain$stim)

prop.test(x = c(1454, 1495), n = c(2787, 2758),alternative = "two.sided")

xtabs (~ Exc$choice + Exc$previousTr)

prop.test(x = c(1454, 1126), n = c(2787, 2217),alternative = "two.sided")

xtabs (~ Sham$choice + Sham$previousTr)

prop.test(x = c(1496, 1260), n = c(2758, 2178),alternative = "two.sided")