######################################################################
############################# user input #############################
###################################################################### 

#---------------------------------------------------------------------
#----------------set directories for data input und data output-------
#---------------------------------------------------------------------

# per default, directories are set for windows machines; directories
# can be adjusted for systems running mac or linux
directoryData = "C:/data.csv" # raw csv data file

# per default, the seperator between values in the rows of data.csv is
# "," - a space/tab, or ";" can be used instead
separator = ","

directoryOutput = "C:/coefficients.csv" # set directory for estimated
                                        # coefficients
directorylogLikBic = "C:/likbic.csv"    # set directory for log-likelihoods

#---------------------------------------------------------------------
#----------------set strategy labels----------------------------------
#---------------------------------------------------------------------

strat1 = "TTB"
strat2 = "EQW"
strat3 = "WADDcorrected"
strat4 = "PCS"
strat5 = "RANDOM"

#---------------------------------------------------------------------
#----------------inactivate/activate code sections--------------------
#---------------------------------------------------------------------

readdata    = 1  # read in data? 1 = yes, 0 = no
mmml        = 1  # do MM-ML? 1 = yes, 0 = no
totalprint  = 1  # print coefficients, log-likelihoods + BICs
                 #  for all participants? 1 = yes, 0 = no
saveoutput  = 1  # save coefficients, log-likelihoods + BICs?
                 # 1 = yes, 0 = no

#---------------------------------------------------------------------
#---output: print selection of coefficients, log-likelihoods and BICs 
#---(by participant and strategy)-------------------------------------
#---------------------------------------------------------------------

participantOUT = 6  # set the number of the participant of interest
strategyOUT    = 4  # set the number of the strategy of interest:
                    # 1 = TTB, 2 = EQW, 3 = WADDcorrected, 4 = PCS,
                    # 5 = RANDOM 

#---------------------------------------------------------------------
#------configure the properties of maximum likelihood algorithm (mle)
#---------------------------------------------------------------------

# set maximum number of iterations for mm-ml; in case the mle
# algorithm does not converge (which terminates program exicution and
# produces an error message), increase the number of iterations

numberOFiterations = 10^7   

# set the criterion of convergence tolerance for the mle algorithm; a
# lower criterion will lead to more computation --> enter
# "help(optim)" in the open R console and hit "ENTER" for a
# description of the mle algorithm; the criterion should not be >
# 10^-9 and not < 0

relconvergencefactor = 10^-15  

########################################################################
#program starts here; beyond this line, no further user input is needed#
########################################################################  

#---------------------------------------------------------------------
#-------------------read in data--------------------------------------
#---------------------------------------------------------------------

if (readdata == 1){

# raw data file is loaded
dat = read.csv(directoryData , sep = separator , header = T)

partic  = dat[,1]  # columns of data.csv are separated into vectors
strat   = dat[,3]  # (see table 3)
type    = dat[,4]
index   = dat[,5]
total   = dat[,6]
miss    = dat[,7]
dv      = dat[,8]

numberOFstrategies   = max(strat)         # the number of strategies
numberOFparticipants = max(partic)        # number of participants
numberOFtasks        = max(total[index])  # number of tasks per type
numberOFtypes        = max(type)          # number of types of tasks

hit1 = dv
# change DV values of 1.1 to 1 
hit1[index==1] = ifelse(dv[index==1]==1.1,1,dv[index==1])
# vector for choices  
choice = matrix(,nrow = max(dat[,1]),ncol =
  numberOFstrategies*numberOFtypes)
# vector for decision times
time = matrix(,nrow = max(dat[,1]),ncol =
  numberOFstrategies*numberOFtasks*numberOFtypes)
# vector for confidence judgments
confidence = matrix(,nrow = max(dat[,1]),ncol =
  numberOFstrategies*numberOFtasks*numberOFtypes)
# vector for the contrast for choices
Tchoice = matrix(,nrow = max(dat[,1]),ncol =
  numberOFstrategies*numberOFtypes)
# vector for the contrast for decsion times
Ttime = matrix(,nrow = max(dat[,1]),ncol =
  numberOFstrategies*numberOFtasks*numberOFtypes)
# vector for the contrast for confidence judgments 
Tconfidence =  matrix(,nrow = max(dat[,1]),ncol =
  numberOFstrategies*numberOFtasks*numberOFtypes)

# loop over participants, strategies, and tasks to store choices
for(i in (1:numberOFparticipants)){
g = 1
for(k in (1:numberOFstrategies)){
for(j in (1:numberOFtypes)){

    for(z in (1:length(dv))){
               
        if(partic[z] == i & strat[z] == k & type[z] == j &
           index[z] == 1 & miss[z] == 0){
            choice[i,g] = hit1[z]
            Tchoice[i,g] = 1
            g = g+1}
        if(partic[z] == i & strat[z] == k & type[z] == j &
           index[z] == 1 & miss[z] == 1){
            choice[i,g] = hit1[z]
            Tchoice[i,g] = 0
            g = g+1}}    
}}

print(paste("choices for participant #",i,"loaded"))
}

# loop over participants, strategies, and tasks to store decision times
for(i in (1:numberOFparticipants)){     
g = 1
for(k in (1:numberOFstrategies)){
for(j in (1:numberOFtypes)){                                  

    for(z in (1:length(dv))){
        if(partic[z] == i & strat[z] == k & type[z] == j &
           index[z] == 2){
            Ttime[i,g] = total[z]
            time[i,g] = hit1[z]
            g = g+1}}
}}

print(paste("decision times for participant #",i,"loaded"))
}

# loop over participants, strategies, and tasks to store confidence judgments
for(i in (1:numberOFparticipants)){    
g = 1
for(k in (1:numberOFstrategies)){
for(j in (1:numberOFtypes)){

    for(z in (1:length(dv))){
        if(partic[z] == i & strat[z] == k & type[z] == j & index[z] == 3){
            Tconfidence[i,g] = total[z]
            confidence[i,g] = hit1[z]
            g = g+1}}
}}

print(paste("confidence judgments for participant #",i,"loaded"))
}

# recode choices from number of strategy inconsistent to consistent choices
choice = abs(choice-numberOFtasks)      
}

#---------------------------------------------------------------------
#-----start values for the coefficients used by the mle algorithm-----
#---------------------------------------------------------------------

# any value > 0 and < 1 is o.k.             
startEPSILON          =  .5                                                             
# a reasonable starting value for MUtime is the observed total mean
# decision times averaged over participants
startMUtime           =  mean(colMeans(time[,1:(numberOFtypes*numberOFtasks)]))         
# a reasonable starting value for SIGMAtime is the observed total
# standard deviation of decision times averaged over participants
startSIGMAtime        =  sd(colMeans(time[,1:(numberOFtypes*numberOFtasks)]))           
# a reasonable starting value for MUtconfidence is the observed total
# mean of confidence judgments averaged over participants
startMUconfidence     =  mean(colMeans(confidence[,1:(numberOFtypes*numberOFtasks)]))   
# a reasonable starting value for SIGMAconfidence is the observed total
# standard deviation of confidence judgments averaged over participants
startSIGMAconfidence  =  sd(colMeans(confidence[,1:(numberOFtypes*numberOFtasks)]))     

# a reasonable starting value for startRtime is the observed standard
# deviation between the means of decision times per type of tasks
# averaged over participants
startRtime = colMeans(time[,1:(numberOFtypes*numberOFtasks)])                           
startRtime = matrix(startRtime,numberOFtasks)
startRtime = sd(colMeans(startRtime))

# a reasonable starting value for startRconfidence is the observed
# standard deviation between the means of confidence judgments per
# type of tasks averaged over participants
startRconfidence = colMeans(confidence[,1:(numberOFtypes*numberOFtasks)])               
startRconfidence = matrix(startRconfidence,numberOFtasks)
startRconfidence = sd(colMeans(startRconfidence))

#--------------------------------------------------------------------
#----------------load the required R library for the mle algorithm---
#--------------------------------------------------------------------

if (mmml == 1){

library(stats4)              

#--------------------------------------------------------------------
#----------------define required output variables--------------------
#--------------------------------------------------------------------

output = numeric()
likbic = numeric()

######################################################################
############ maximum likelihood estimation (mle) starts here #########
###################################################################### 

# loop over all participants
for (participants in 1 : numberOFparticipants){     

# general comment: for the strategies discussed in the paper, there
# are 7 (TTB and PCS), 6 (EQW and WADDcorrected), and 4 (RANDOM)
# coefficients that need to be estimated by the mle algorithm. We
# therefore need to specifiy 3 (= 7/6/4 coefficients) * 3 (= choices,
# decision times, and confidence judgments) functions in total. In the
# following, the program is commented for the first strategies (TTB
# and PCS) and only notable differences are commented for the
# remaining strategies (= EQW, WADDcorrected, and RANDOM).

#-----------------------------------------------------------------------------
#------------mle for the strategies 'Take The Best' (= TTB) + ----------------
#-------------'Parallel Constraint Satisfaction' (= PCS)----------------------
#-----------------------------------------------------------------------------

# a and b modify the epsilons for choices. For some combinations of
# strategies and type of tasks, e. g., EQW and type of tasks 2, both
# options (A,B) are in line with the predicted choice. For these
# combinations epsilon = .5. For the implementation, we define
# (1-epsilon) = (1-epsilon)*a+b.  For combinations where a choice of A
# or B is in line with the strategy a = 0 and b = .5, for all other
# cases a = 1 and b = 0.

a = Tchoice;
b = ifelse(Tchoice == 1,0,1/2)  

# number of tasks per type are stored in a vector
numberoftaskspertype = rep(numberOFtasks,numberOFstrategies * numberOFtypes)      

for (strategy in 1 : 2){      # loop over strategy 1 and 2 

# indexes are set for each strategy

if (strategy == 1){

i = 1;
j = NCOL(choice)/numberOFstrategies;
k = 1;
l = (NCOL(time)/numberOFstrategies); 
}

if (strategy == 2){

i = 3*(NCOL(choice)/numberOFstrategies)+1;
j =  4*(NCOL(choice)/numberOFstrategies);
k = 3*(NCOL(time)/numberOFstrategies)+1; 
l = 4*(NCOL(time)/numberOFstrategies);
}

# It is more convenient to maximize log-likelihood functions than
# likelihood functions (likelihoods tend to be small).

# function for the log-likelihoods for choices (see the first part of
# equation 1)
Choice = function(epsilon = startEPSILON)      
-sum(dbinom(choice[participants,i:j], numberoftaskspertype[i:j],
            prob = ((1-epsilon)*a[participants,i:j]+b[participants,i:j]),
            log = TRUE));
fit1 = mle(Choice, method = "BFGS",control = list(maxit = numberOFiterations,
            reltol = relconvergencefactor));

# function for the log-likelihoods for decision times (see the second
# part of equation 1)
Time = function(mueTIME = startMUtime, sigmaTIME = startSIGMAtime,
  RTIME = startRtime)-sum(dnorm(time[participants,k:l],
    mean = (mueTIME+(Ttime[participants,k:l]*RTIME)),
    sd = sigmaTIME,log = TRUE));
fit2 = mle(Time, method = "BFGS",
  control = list(maxit = numberOFiterations,reltol = relconvergencefactor));
 
# function for the log-likelihoods for confidence judgments (see the
# third part of equation 1)
Confidence = function(mueCONF = startMUconfidence,
  sigmaCONF = startSIGMAconfidence,RCONF = startRconfidence) 
-sum(dnorm(confidence[participants,k:l],
           mean = (mueCONF+(Tconfidence[participants,k:l]*RCONF)),
           sd = sigmaCONF,log = TRUE));    
fit3 = mle(Confidence, method= "BFGS",
  control = list(maxit = numberOFiterations,reltol = relconvergencefactor));

# sum up log-likelihoods for choices, decision times, and confidence judgments
y = (logLik(fit1)+logLik(fit2)+logLik(fit3))[1]   

# calculate the corresponding BIC values
BICs = (-2*(y)+log(3*numberOFtypes)*7)[1]         

# the following variables (v1 to v28) store the coefficients and test
# statistics: 'epsilon' = v1,'z epsilon' = v2,'p epsilon' = v3,
# 'epsilonci' = v4, 'mueTIME' = v5,'z mueTime' = v6,'p mueTime' = v7,
# 'mueTIMEci' = v8,'sigmaTIME' = v9,'z sigmaTIME' = v10,'p
# sigmaTIME' = v11,'sigmaTIMEci' = v12, 'RTIME'=v13,'z RTIME' = v14,
# 'p RTIME' = v15,'RTIMEci' = v16,'mueCONF' = v17,'z mueCONF' = v18,
# 'p mueCONF' = v19,'mueCONFci' = v20, 'sigmaCONF' = v21,'z sigmaCONF' =
# v22,'p sigmaCONF' = v23,'sigmaCONFci' = v24, 'RCONF' = v25,'z RCONF'
# = v26,'p RCONF' = v27, and 'RCONFci' = v28

v1 =  coef(summary(fit1))                              
v2 =  v1[1]/v1[2]   
v3 =  2*(1-pnorm(abs(v2)))
v4 =  c(v1[1]-v1[2]*1.96,v1[1]+v1[2]*1.96)
v5 =  coef(summary(fit2))[1,]
v6 =  v5[1]/v5[2]
v7 =  2*(1-pnorm(abs(v6)))
v8 =  confint(fit2)[1,]
v9 =  coef(summary(fit2))[2,]
v10 = v9[1]/v9[2]
v11 = 2*(1-pnorm(abs(v10)))
v12 = confint(fit2)[2,]
v13 = coef(summary(fit2))[3,]
v14 = v13[1]/v13[2]
v15 = 2*(1-pnorm(abs(v14)))
v16 = confint(fit2)[3,]
v17 = coef(summary(fit3))[1,]
v18 = v17[1]/v17[2]
v19 = 2*(1-pnorm(abs(v18))) 
v20 = confint(fit3)[1,]
v21 = coef(summary(fit3))[2,]
v22 = v21[1]/v21[2]
v23 = 2*(1-pnorm(abs(v22)))
v24 = confint(fit3)[2,]
v25 = coef(summary(fit3))[3,]
v26 = v25[1]/v25[2]
v27 = 2*(1-pnorm(abs(v26)))
v28 = confint(fit3)[3,]

# results are stored in the list TTB
if (strategy == 1){
  TTB = list(y,'BIC' = BICs,'epsilon' = v1,'z' = v2,'p' =
  v3,'epsilonci' = v4,'mueTIME' = v5,'z' = v6,'p' = v7,'mueTIMEci' =
  v8,'sigmaTIME' = v9,'z' = v10,'p' = v11,'sigmaTIMEci' =
  v12,'RTIME'=v13,'z' = v14,'p' = v15,'RTIMEci' = v16,'mueCONF' =
  v17,'z' = v18,'p' = v19,'mueCONFci' = v20,'sigmaCONF' = v21,'z' =
  v22,'p' = v23,'sigmaCONFci' = v24,'RCONF' = v25,'z' = v26,'p' =
  v27,'RCONFci' = v28)
}

# results are stored in the list PCS
if (strategy == 2){
  PCS = list(y,'BIC' = BICs,'epsilon' = v1,'z' = v2,'p' =
  v3,'epsilonci' = v4,'mueTIME' = v5,'z' = v6,'p' = v7,'mueTIMEci' =
  v8,'sigmaTIME' = v9,'z' = v10,'p' = v11,'sigmaTIMEci' = v12,'RTIME'
  = v13,'z' = v14,'p' = v15,'RTIMEci' = v16,'mueCONF' = v17,'z' =
  v18,'p' = v19,'mueCONFci' = v20,'sigmaCONF' = v21,'z' = v22,'p' =
  v23,'sigmaCONFci' = v24,'RCONF' = v25,'z' = v26,'p' = v27,'RCONFci'
  = v28)
}};

#-----------------------------------------------------------------------------
#-------mle for the strategies 'Equal Weight' (= EQW) + ----------------------
# -------'Weighted Additive Corrected' (WADDcorrected)------------------------
#-----------------------------------------------------------------------------

# only 6 coefficients need to be estimated

for (strategy in 3 : 4){

if (strategy == 3){

i = 1*(NCOL(choice)/numberOFstrategies)+1;
j =  2*(NCOL(choice)/numberOFstrategies);
k = 1*(NCOL(time)/numberOFstrategies)+1; 
l = 2*(NCOL(time)/numberOFstrategies); 

}

if (strategy == 4){

i = 2*(NCOL(choice)/numberOFstrategies)+1;
j = 3*(NCOL(choice)/numberOFstrategies);
k = 2*(NCOL(time)/numberOFstrategies)+1;
l = 3*(NCOL(time)/numberOFstrategies);

}

Choice = function(epsilon = startEPSILON)-
  sum(dbinom(choice[participants,i:j], numberoftaskspertype[i:j],
             prob = ((1-epsilon)*a[participants,i:j]+b[participants,i:j]),
             log = TRUE));
fit1 = mle(Choice, method = "BFGS",
  control = list(maxit = numberOFiterations,reltol = relconvergencefactor));

# there is no RTIME that needs to be estimated for EQW and WADDcorrected
Time = function(mueTIME = startMUtime, sigmaTIME = startSIGMAtime)      
-sum(dnorm(time[participants,k:l], mean = (mueTIME),sd = sigmaTIME,log = TRUE));
fit2 = mle(Time, method = "BFGS",control = list(maxit = numberOFiterations,
             reltol = relconvergencefactor));

Confidence = function(mueCONF = startMUconfidence,
  sigmaCONF = startSIGMAconfidence,RCONF = startRconfidence)
-sum(dnorm(confidence[participants,k:l],
           mean = (mueCONF+(Tconfidence[participants,k:l]*RCONF)),
           sd = sigmaCONF,log = TRUE));
fit3 = mle(Confidence, method = "BFGS",
  control = list(maxit = numberOFiterations,reltol = relconvergencefactor));

y = (logLik(fit1)+logLik(fit2)+logLik(fit3))[1]
BICs = (-2*(y)+log(3*numberOFtypes)*6)[1]

v1 =  coef(summary(fit1))
v2 =  v1[1]/v1[2]   
v3 =  2*(1-pnorm(abs(v2)))
v4 =  c(v1[1]-v1[2]*1.96,v1[1]+v1[2]*1.96)
v5 =  coef(summary(fit2))[1,]
v6 =  v5[1]/v5[2]
v7 =  2*(1-pnorm(abs(v6)))
v8 =  confint(fit2)[1,]
v9 =  coef(summary(fit2))[2,]
v10 = v9[1]/v9[2]
v11 = 2*(1-pnorm(abs(v10)))
v12 =  confint(fit2)[2,]
v13 = coef(summary(fit3))[1,]
v14 = v13[1]/v13[2]
v15 = 2*(1-pnorm(abs(v14))) 
v16 = confint(fit3)[1,]
v17 = coef(summary(fit3))[2,]
v18 = v17[1]/v17[2]
v19 = 2*(1-pnorm(abs(v22)))
v20 = confint(fit3)[2,]
v21 = coef(summary(fit3))[3,]
v22 = v21[1]/v21[2]
v23 = 2*(1-pnorm(abs(v22)))
v24 = confint(fit3)[3,]

if (strategy == 3){
# results are stored in the list EQW 
EQW = list(y,'BIC' =
BICs,'epsilon' = v1,'z' = v2,'p' = v3,'epsilonci' = v4,'mueTIME' =
v5,'z' = v6,'p' = v7,'mueTIMEci' = v8,'sigmaTIME' = v9,'z' = v10,'p' =
v11,'sigmaTIMEci' = v12,'mueCONF' = v13,'z' = v14,'p' =
v15,'mueCONFci' = v16,'sigmaCONF' = v17,'z' = v18,'p' =
v19,'sigmaCONFci' = v20,'RCONF' = v21,'z' = v22,'p' = v23,'RCONFci' =
v24)
}

if (strategy == 4){
# results are stored in the list WADDcorrected 
WADDcorrected =
list(y,'BIC' = BICs,'epsilon' = v1,'z' = v2,'p' = v3,'epsilonci' =
v4,'mueTIME' = v5,'z' = v6,'p' = v7,'mueTIMEci' = v8,'sigmaTIME' =
v9,'z' = v10,'p' = v11,'sigmaTIMEci' = v12,'mueCONF' = v13,'z' =
v14,'p' = v15,'mueCONFci' = v16,'sigmaCONF' = v17,'z' = v18,'p' =
v19,'sigmaCONFci' = v20,'RCONF' = v21,'z' = v22,'p' = v23,'RCONFci' =
v24)
}}

#-----------------------------------------------------------------------------
#-------------mle for the strategy 'RANDOM'-----------------------------------
#-----------------------------------------------------------------------------

# only 4 coefficients need to be estimated

i = 1;
j = NCOL(choice)/numberOFstrategies;
k = 1;
l = (NCOL(time)/numberOFstrategies); 

Choice = function(epsilon = 0.5) -sum(dbinom(choice[participants,i:j],
numberoftaskspertype[i:j], prob = 1-epsilon, log = TRUE)); fit1 =
mle(Choice, method = "BFGS",control = list(maxit =
numberOFiterations,reltol = relconvergencefactor),fixed = list(epsilon
= 0.5)); #epsilon = .5 for RANDOM  

# there is no RTIME that needs to be estimated for RANDOM
Time = function(mueTIME = startMUtime, sigmaTIME = startSIGMAtime)      
-sum(dnorm(time[participants,k:l], mean = mueTIME,sd = sigmaTIME,log = TRUE));
fit2 = mle(Time, method = "BFGS",control = list(maxit = numberOFiterations,
                                   reltol = relconvergencefactor));
                                   
# there is no RCONF that needs to be estimated for RANDOM 
Confidence =
function(mueCONF = startMUconfidence, sigmaCONF =
startSIGMAconfidence) -sum(dnorm(confidence[participants,k:l], mean =
mueCONF,sd = sigmaCONF,log = TRUE)); fit3 = mle(Confidence, method =
"BFGS",control = list(maxit = numberOFiterations,reltol =
relconvergencefactor));

RANDOMy = (logLik(fit1)+logLik(fit2)+logLik(fit3))[1]
RANDOMBIC = (-2*(RANDOMy)+log(3*numberOFtypes)*4)[1]

v1 =  coef(summary(fit2))[1,]
v2 =  v1[1]/v1[2]
v3 =  2*(1-pnorm(abs(v2)))
v4 =  confint(fit2)[1,]
v5 =  coef(summary(fit2))[2,]
v6 = v5[1]/v5[2]
v7 = 2*(1-pnorm(abs(v6)))
v8 =  confint(fit2)[2,]
v9 = coef(summary(fit3))[1,]
v10 = v9[1]/v9[2]
v11 = 2*(1-pnorm(abs(v10))) 
v12 = confint(fit3)[1,]
v13 = coef(summary(fit3))[2,]
v14 = v13[1]/v13[2]
v15 = 2*(1-pnorm(abs(v14)))
v16 = confint(fit3)[2,]

# results are stored in the list RANDOM
RANDOM = list(RANDOMy,'BIC' = RANDOMBIC,'mueTIME' = v1,'z' = v2,'p' =
v3,'mueTIMEci' = v4,'sigmaTIME' = v5,'z' = v6,'p' = v7,'sigmaTIMEci' =
v8,'mueCONF' = v9,'z' = v10,'p' = v11,'mueCONFci' = v12,'sigmaCONF' =
v13,'z' = v14,'p' = v15,'sigmaCONFci' = v16)

######################################################################
############################# Output #################################
###################################################################### 

# for each strategy, the original list (TTB, EQW, ...) is unpacked and
# stored in a more convenient matrix representation
TTBout = unlist(TTB)      
TTBout = t(matrix((TTBout[3:length(TTBout)]),6,7))

EQWout = unlist(EQW)
EQWout = t(matrix((EQWout[3:length(EQWout)]),6,6))

WADDcorrectedout = unlist(WADDcorrected)
WADDcorrectedout = t(matrix((WADDcorrectedout[3:length(WADDcorrectedout)]),6,6))

PCSout = unlist(PCS)
PCSout = t(matrix((PCSout[3:length(PCSout)]),6,7))

RANDOMout = unlist(RANDOM)
RANDOMout = t(matrix((RANDOMout[3:length(RANDOMout)]),6,4))

# labels are set for the output
rnames1 = c("epsilon","mu_Time","sigma_Time","R_Time","mu_Conf","sigma_Conf",
  "R_Conf")
rnames2 = c("epsilon","mu_Time","sigma_Time","mu_Conf","sigma_Conf","R_Conf")
rnames3 = c("mu_Time","sigma_Time","mu_Conf","sigma_Conf")

row.names(TTBout) = rnames1
row.names(EQWout) = rnames2
row.names(WADDcorrectedout) = rnames2
row.names(PCSout) = rnames1
row.names(RANDOMout)=rnames3

# the final output matrix is set
output = rbind(output,cbind(rep(participants,26),c(rep(1,NROW(TTBout)),
  rep(2,NROW(EQWout)),rep(3,NROW(WADDcorrectedout)),rep(4,NROW(PCSout)),
  rep(5,NROW(RANDOMout))),rbind(TTBout,EQWout,WADDcorrectedout,PCSout,RANDOMout))) 
colnames(output) = c("participant","strategy","Coef.","Std. Err.","z",
          "P>|z|","CI 2.5%","CI 97.5%")

likbic = rbind(likbic,cbind(rep(participants,2),c(1,2),cbind(unlist(TTB)[1:2],
  unlist(EQW)[1:2],unlist(WADDcorrected)[1:2],unlist(PCS)[1:2]),c(RANDOMy,RANDOMBIC)))
colnames(likbic) = c("participant","Lik/Bic",strat1,strat2,strat3,strat4,strat5)

print(paste("MM-ML for participant #", participants, "done") )
}

rownames(likbic) = rep(c("logLik","BIC"),participants)
}

#-----------------------------------------------------------------------------
#------------------print output-----------------------------------------------
#-----------------------------------------------------------------------------

if (totalprint == 1){

print(output)     # coefficients for all participants + strategies
print(likbic)     # log-likelihoods and BICs for all participants + strategies

likbic[likbic[,2] == 1,]   # log-likelihoods for all participants + strategies
likbic[likbic[,2] == 2,]   # BICs for all participants + strategies
}

#-----------------------------------------------------------------------------
#----------print a subset of the output (set under user input)----------------
#-----------------------------------------------------------------------------
                                   
likbic[likbic[,1] == participantOUT,]
output[output[,1] == participantOUT & output[,2] == strategyOUT,]

#-----------------------------------------------------------------------------
#-------------------------save output (set under user input)------------------
#-----------------------------------------------------------------------------

if (saveoutput == 1){
write.csv(output,file = directoryOutput,row.names = FALSE)
write.csv(likbic,file = directorylogLikBic,row.names = FALSE)
}