### function name, variables, and defaults MMML = function(partic = 1:max(dat$PARTICIPANT),stratname="NONAME",expectedchoice = "none",contrasttime = "none" ,contrastconfidence = "none", saveoutput = 1,directoryData = "C:/datafunction.csv", directoryOutput = "C:/output.csv",separator = ",",numberOFiterations = 10^9, relconvergencefactor = 10^-15){ ### load required library for maximum likelihood estimation (mle) library(stats4) ### set control values for mle controlParameterMLE = list(maxit = numberOFiterations,reltol = relconvergencefactor) ### read in data dat = read.csv(directoryData,sep = separator,header = T) ### define variables and set default values guessing = 0 # guessing for choices zerocontrasttime = 0 # zero contrast for decision times zerocontrastconfidence = 0 # zero contrast for confidence judgments counterGUESSING = 0 # needed to check for guessing memory = numeric() # output is stored in the variable memory choice = numeric() ### check if there is only guessing for choices if(expectedchoice[1] != "none"){ for (i in 1 : length(expectedchoice)){ if(expectedchoice[i] == 1/2){ counterGUESSING = counterGUESSING + 1 } } if (counterGUESSING==length(expectedchoice)){ guessing = 1 } } ### check if there is a zero contrast for decision times if (contrasttime[1]!="none"){ if (sum(abs(contrasttime)) == 0 | sd(contrasttime)==0){ zerocontrasttime = 1 contrasttime = rep(0,length(contrasttime)) } } ### check if there is a zero contrast for confidence judgments if (contrastconfidence[1]!="none"){ if (sum(abs(contrastconfidence)) == 0 | sd(contrastconfidence)==0){ zerocontrastconfidence = 1 contrastconfidence = rep(0,length(contrastconfidence)) } } ### determine the number of parameters that need to be ### estimated (needed to calculate BICs) numberOFfreeParameters = 7 # full model if(guessing == 1){ numberOFfreeParameters=numberOFfreeParameters-1 } if(expectedchoice[1] == "none"){ numberOFfreeParameters=numberOFfreeParameters-1 } if(contrasttime[1] == "none"){ numberOFfreeParameters=numberOFfreeParameters-3 } if(contrastconfidence[1] == "none"){ numberOFfreeParameters=numberOFfreeParameters-3 } if(zerocontrasttime == 1){ numberOFfreeParameters=numberOFfreeParameters-1 } if(zerocontrastconfidence == 1){ numberOFfreeParameters=numberOFfreeParameters-1 } ### change contrasts for decision times and confidence ### judgments --> range of contrasts = 1 and sum of contrasts = 0 if(contrasttime[1] != "none" & zerocontrasttime !=1){ contrasttime = contrasttime-(sum(contrasttime)/length(contrasttime)) rangetime=range(contrasttime)[2]-range(contrasttime)[1] contrasttime = contrasttime/rangetime } if(contrastconfidence[1] != "none" & zerocontrastconfidence != 1){ contrastconfidence=contrastconfidence- (sum(contrastconfidence)/length(contrastconfidence)) rangeconfidence=range(contrastconfidence)[2]- range(contrastconfidence)[1] contrastconfidence = contrastconfidence/rangeconfidence } ### loop participants for(numbpartic in 1 : length(partic)){ ### define variables numberOFtypes = max(dat$TYPE[dat$PARTICIPANT==partic[numbpartic]]) numberOFtasksPerType = numeric() for(loopnumberOFtasksPerType in 1 : numberOFtypes){ numberOFtasksPerType[loopnumberOFtasksPerType] = NROW(dat[dat$TYPE == loopnumberOFtasksPerType & dat$PARTICIPANT ==partic[numbpartic],]) } logliks = c(0,0,0) numberOFmeasures = 0 stratOutput = numeric() epsilon = numeric() messageerror = "" errortermneeded = 1 # see *** below --> is it necessary # to add/subtract an error term? ### prepare data vector for choices if(expectedchoice[1] != "none"){ for(looptype in 1 : numberOFtypes){ choice[looptype] = sum(dat$ACHOICES[dat$PARTICIPANT == partic[numbpartic] & dat$TYPE == looptype]) } OnEpsilon = ifelse(expectedchoice==1/2,0,1) # see the function for choices below --> if # the strategy predicts option A or B, epsilon # needs to be fixed at 0.5 FixedEpsilon = ifelse(expectedchoice==1/2,1/2,0) ### If all choices are (resp. are not) in line # *** ### with the expected choices, it is necessary ### for mle convergence to add ### (resp. subtract) 1 strategy inconsistent ### (resp. consistent) choice. It is further ### necessary that the 1 strategy inconsistent ### (resp. consistent) choice is added to ### (resp. subtracted from) a type of tasks ### that does NOT predict guessing (i.e, fixed ### epsilon = .5). for (looperror in 1 : length(numberOFtasksPerType)){ if(choice[looperror]!=numberOFtasksPerType[looperror] & choice[looperror]!=0 & expectedchoice[looperror] != 1/2){ errortermneeded = 0 } } if (guessing != 1 & errortermneeded == 1){ typeerror = which(expectedchoice!=1/2)[1] if (choice[typeerror] == 0){ choice[typeerror] = 1 messageerror = paste("ERROR TERM: 1 A CHOICE ADDED TO TYPE OF TASKS",typeerror) }else{ choice[typeerror] = numberOFtasksPerType[typeerror]-1 messageerror = paste("ERROR TERM: 1 A CHOICE SUBTRACTED FROM TYPE OF TASKS", typeerror) } } ### if option B is predicted, choices and expected choices need ### to be inverted expectedchoiceINVERTED = ifelse(expectedchoice == 1/2, numberOFtasksPerType,expectedchoice) indexItemsInverted = which(expectedchoiceINVERTED == 0) for (loopinversion in 1 : length(indexItemsInverted)){ choice[indexItemsInverted[loopinversion]] = numberOFtasksPerType[indexItemsInverted[loopinversion]] - choice[indexItemsInverted[loopinversion]] expectedchoiceINVERTED[indexItemsInverted[loopinversion]] = numberOFtasksPerType[indexItemsInverted[loopinversion]] } numberOFmeasures = numberOFmeasures + 1 } ### prepare data vector for decision times if(contrasttime[1] != "none"){ time = dat$DECTIMES[dat$PARTICIPANT == partic[numbpartic]] typesTime = dat$TYPE[dat$PARTICIPANT == partic[numbpartic]] Ttime = contrasttime[typesTime] numberOFmeasures = numberOFmeasures + 1 } ### prepare data vector for confidence judgments if(contrastconfidence[1] != "none"){ confidence = dat$CONFJUDGMENTS[dat$PARTICIPANT == partic[numbpartic]] typesConfidence = dat$TYPE[dat$PARTICIPANT == partic[numbpartic]] Tconfidence = contrastconfidence[typesConfidence] numberOFmeasures = numberOFmeasures + 1 } ################################################################### ############### define density functions ########################## ################################################################### ### choices if(expectedchoice[1] != "none"){ if(guessing!=1){ Choice = function(epsil){-sum(dbinom(choice, expectedchoiceINVERTED, prob = (((1-epsil) * OnEpsilon) + FixedEpsilon) , log = TRUE))}; }else{ Choice = function(epsil){-sum(dbinom(choice, expectedchoiceINVERTED, prob = 1-epsil , log = TRUE))}; # function for the RANDOM model } } ### decision times if(zerocontrasttime != 1){ Time = function(mu_Time = 8,sigma_Time, R_Time){-sum(dnorm(time, mean = (mu_Time+(Ttime*R_Time)), sd = sigma_Time,log = TRUE))}; }else{ Time = function(mu_Time, sigma_Time){-sum(dnorm(time, mean = mu_Time,sd = sigma_Time,log = TRUE))}; } ### confidence judgments if(zerocontrastconfidence != 1){ Confidence = function(mu_Conf, sigma_Conf,R_Conf){-sum(dnorm(confidence, mean = (mu_Conf +(Tconfidence*R_Conf)),sd = sigma_Conf,log = TRUE))}; }else{ Confidence = function(mu_Conf, sigma_Conf){-sum(dnorm(confidence, mean = mu_Conf ,sd = sigma_Conf,log = TRUE))}; } ################################################################### ############## do mle ############################################# ################################################################### ### set starting values for mle start_epsilon = 0.5 if(contrasttime[1] != "none"){ start_mu_Time = mean(time) start_sigma_Time = sd(time) start_R_Time = sd(tapply(time,typesTime,mean)) } if(contrastconfidence[1] != "none"){ start_mu_Conf = mean(confidence) start_sigma_Conf = sd(confidence) start_R_Conf = sd(tapply(confidence,typesConfidence,mean)) } ### mle for choices if(expectedchoice[1] != "none"){ if(guessing != 1){ fit1 = mle(Choice, start=list(epsil= start_epsilon), method = "BFGS",control = substitute(controlParameterMLE)); }else{ fit1 = mle(Choice, start=list(epsil= start_epsilon), method = "BFGS",control = substitute(controlParameterMLE), fixed = list(epsil = 0.5)); } logliks[1] = logLik(fit1) } ### mle for decision times if(contrasttime[1] != "none"){ if(zerocontrasttime != 1){ startTIME = list(mu_Time = start_mu_Time, sigma_Time = start_sigma_Time,R_Time = start_R_Time) }else{ startTIME = list(mu_Time = start_mu_Time, sigma_Time = start_sigma_Time) } fit2 = mle(Time, start = startTIME, method ="BFGS", control = substitute(controlParameterMLE)); logliks[2] = logLik(fit2) } ### mle for confidence judgments if(contrastconfidence[1] != "none"){ if(zerocontrastconfidence != 1){ startCONFIDENCE = list(mu_Conf=start_mu_Conf, sigma_Conf=start_sigma_Conf,R_Conf=start_R_Conf) }else{ startCONFIDENCE = list(mu_Conf=start_mu_Conf, sigma_Conf=start_sigma_Conf) } fit3 = mle(Confidence,start=startCONFIDENCE, method ="BFGS", control =substitute(controlParameterMLE)); logliks[3] = logLik(fit3)[1] } ################################################################### ######### calculate log-liks, BICs, coefficients, z-scores, ####### ######### significance levels, and confidence intervals ########## ################################################################### ### log-likelihoods logliks = sum(logliks) ### BICs BICs = (-2*(logliks)+log(numberOFmeasures*numberOFtypes)* numberOFfreeParameters)[1] ### choices if(expectedchoice[1] != "none" & guessing != 1){ epsilon = coef(summary(fit1)) epsilon = c(epsilon ,coef(summary(fit1))[,1]/coef(summary(fit1))[,2]) epsilon = c(epsilon,2*(1-pnorm(abs(epsilon[3])))) epsilon = c(epsilon, epsilon[1]-epsilon[2]*1.96, epsilon[1]+epsilon[2]*1.96) ### set labels if MM-ML is only based on choices if(contrasttime[1] == "none" & contrastconfidence[1] == "none"){ names(epsilon) = c("epsilon","Std. Error epsilon", "z epsilon","P>|z| epsilon","2.5 % CI epsilon", "97.5 % CI epsilon") } } ### decision times if(contrasttime[1] != "none"){ outputSTRATtime = coef(summary(fit2)) outputSTRATtime = cbind(outputSTRATtime, "z"=coef(summary(fit2))[,1]/ coef(summary(fit2))[,2]) outputSTRATtime = cbind(outputSTRATtime, "P>|z|"=2*(1-pnorm(abs(outputSTRATtime[,3])))) outputSTRATtime = cbind(outputSTRATtime,confint(fit2)) stratOutput = outputSTRATtime } ### confidence judgments if(contrastconfidence[1] != "none"){ outputSTRATconfidence = coef(summary(fit3)) outputSTRATconfidence = cbind(outputSTRATconfidence, "z"=coef(summary(fit3))[,1]/coef(summary(fit3))[,2]) outputSTRATconfidence = cbind(outputSTRATconfidence, "P>|z|"=2*(1-pnorm(abs(outputSTRATconfidence[,3])))) outputSTRATconfidence = cbind(outputSTRATconfidence,confint(fit3)) stratOutput = rbind(stratOutput,outputSTRATconfidence) } stratOutput = rbind(epsilon,stratOutput) ################################################################### ################ Output ########################################### ################################################################### print("###################################################################") print (paste("|","participant =",partic[numbpartic],"| strategy =",stratname,"|")) print (paste("| Log-Likelihood =",round(logliks,2),"| BIC =",round(BICs,2),"|")) if(messageerror!=""){ print(messageerror) } print(stratOutput) print("###################################################################") if (is.matrix(stratOutput)==TRUE){ memory = rbind(memory,cbind("participant"=partic[numbpartic], "strategy name" = stratname,"log-lik"=logliks,"BIC"=BICs,stratOutput)) }else{ memory = rbind(memory,c("participant"=partic[numbpartic], "strategy name" = stratname,"log-lik"=logliks,"BIC"=BICs,stratOutput)) } } # end participant loop ### save data if(saveoutput == 1){ write.csv(memory , file = directoryOutput,row.names = TRUE) } } # end function