Diffusion Models of Pragmatic Inference Tasks
Dan Yurovsky & Erica Yoon
January 07, 2014
Data Loading
Load required Libraries
rm(list=ls())
library(ggplot2)
library(dplyr)
library(tidyr)
library(magrittr)
library(directlabels)
library(RWiener)Fixing up the wiener_plot function to deal with missing data
wiener.plot <- function (dat)
{
rt = as.double(dat$q)
rc = as.numeric(as.factor(dat$resp))
dpos = tryCatch(density(rt[rc == 1], from = 0),error=function(e) NA)
dneg = tryCatch(density(rt[rc == 2], from = 0),error=function(e) NA)
maxt = max(pretty(max(rt)))
maxd <- NA
if(is.na(dpos[1])){
maxd <- max(dneg$y)
} else if(is.na(dneg[1])){
maxd <- max(dpos$y)
} else {
maxd <- max(dpos$y, dneg$y)
}
par(mar = c(0, 5, 0, 0), mfcol = c(2, 1), ask = FALSE)
plot(dpos, xlim = c(0, maxt), ylim = c(0, maxd), las = 2,
lwd = 2, col = "green3", main = "", ylab = "", ask = FALSE)
rug(rt[rc == 1], col = "green3")
mtext("Density of positive responses", side = 2, line = 4,
cex = 0.8)
plot(dneg, xlim = c(0, maxt), ylim = c(maxd, 0), las = 2,
lwd = 2, col = "red", main = "", ylab = "", ask = FALSE)
mtext("Density of negative responses", side = 2, line = 4,
cex = 0.8)
rug(rt[rc == 2], col = "red", side = 3)
}Read in adult data
data <- read.csv('150105-simpimpsc-turk_complete.csv') %>%
filter(trialtype %in% c("control_double", "inference")) %>%
group_by(WorkerId,trialnum) %>%
mutate(side = if((trialtype == "control_double") & (correct == "Y")) "more"
else if((trialtype == "control_double") & (correct == "N")) "less"
else if((trialtype == "inference") & (correct == "N")) "more"
else "less") %>%
mutate(resp = as.character(factor(side,levels=c("more","less"),
labels=c("upper","lower")))) %>%
mutate(q = rt/1000) %>%
group_by() %>%
select(q,trialtype,itemnum,resp) %>%
# filter rt, remove rt outside of stdev*2
filter(log(q) < mean(log(q)) + 2 * sd(log(q)),
log(q) > mean(log(q)) - 2 * sd(log(q)))
control.2v1.data <- filter(data,trialtype == "control_double",
itemnum=="2vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(control.2v1.data)
control.3v1.data <- filter(data,trialtype == "control_double",
itemnum=="3vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(control.3v1.data)
inference.2v1.data <- filter(data,trialtype == "inference",
itemnum=="2vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(inference.2v1.data)
inference.3v1.data <- filter(data,trialtype == "inference",
itemnum=="3vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(inference.3v1.data)
Compute parameters for each condition separately
optim.control.2v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=control.2v1.data, method="Nelder-Mead")
optim.control.3v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=control.3v1.data, method="Nelder-Mead")
optim.inference.2v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=inference.2v1.data, method="Nelder-Mead")
optim.inference.3v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=inference.3v1.data, method="Nelder-Mead")
indiv.pars <- as.data.frame(rbind(optim.control.2v1$par,
optim.control.3v1$par,
optim.inference.2v1$par,
optim.inference.3v1$par))
names(indiv.pars) <- c("Separation","Non.Decision","Bias","Drift")
indiv.pars$Condition <- c("Control", "Control", "Inference", "Inference")
indiv.pars$Item.Num<- c(2, 3, 2, 3)
print(indiv.pars)## Separation Non.Decision Bias Drift Condition Item.Num
## 1 1.786754 0.4291568 0.4340421 1.1447030 Control 2
## 2 2.996473 0.2842518 0.3047898 2.2676962 Control 3
## 3 1.974870 0.4363610 0.5747067 -0.8593789 Inference 2
## 4 2.632747 0.3087888 0.5847111 -1.4294429 Inference 3Pin parameters for item.num together
# Function for tying parameters together across conditions
many_drifts <- function(x, datlist) {
l = 0
for (c in 1:length(datlist)) {
l = l + wiener_deviance(x[c(1, 2, 3, c+3)], datlist[[c]])
}
return(l)
}
datlist.2v1 <- list(control.2v1.data, inference.2v1.data)
datlist.3v1 <- list(control.3v1.data, inference.3v1.data)
# use nlm to estimate parameters
optim.2v1 <- optim(p=c(1, .1, .1, 1, 1), many_drifts,
dat=datlist.2v1,method="Nelder-Mead")
optim.3v1 <- optim(p=c(1, .1, .1, 1, 1), many_drifts,
dat=datlist.3v1,method="Nelder-Mead")
joint.pars <- as.data.frame(rbind(optim.2v1$par,optim.3v1$par))
names(joint.pars) <- c("Separation","Non.Decision","Bias",
"Control", "Inference")
joint.pars <- gather(joint.pars,Condition,Drift,Control:Inference)
joint.pars$Item.Num <- c(2,3,2,3)
joint.pars %<>% select(Separation,Non.Decision,Bias,
Drift,Condition,Item.Num) %>%
arrange(Condition,Item.Num)
print(joint.pars)## Separation Non.Decision Bias Drift Condition Item.Num
## 1 1.861372 0.4442728 0.5132973 0.9776558 Control 2
## 2 2.706737 0.3790312 0.4433851 1.8211632 Control 3
## 3 1.861372 0.4442728 0.5132973 -0.6758614 Inference 2
## 4 2.706737 0.3790312 0.4433851 -1.1623999 Inference 3Pin non-decision time for a group, pin boundaries and bias within a condition
# Function for tying parameters together across conditions
multi.loss.func <- function(x, datmat) {
l = 0
num.conds <- length(datmat)
num.types <- length(datmat[[1]])
# model takes (sep, non-decision, bias, drift)
# x has (non-decision, sep-1, bias-1, drift-11, drift-12,
# sep-2, bias-2, drift-21, drift-22)
for (cond in 1:num.conds) {
for(type in 1:num.types){
l = l + wiener_deviance(x[c(2+(cond-1)*(2+num.types), 1,
3+(cond-1)*(2+num.types),
3+(cond-1)*(2+num.types)+type)],
datmat[[cond]][[type]])
}
}
return(l)
}
#Function for reconstructing the parameter matrix
make.par.mat <- function(optim.output) {
pars <- as.data.frame(t(optim.output$par))
names(pars) <- c("NonDecision","Separation.2v1","Bias.2v1",
"Drift.Control.2v1","Drift.Inference.2v1",
"Separation.3v1","Bias.3v1",
"Drift.Control.3v1","Drift.Inference.3v1")
par.mat <- expand.grid(Condition = c("2v1", "3v1"),
Trial.Type = c("Control", "Inference"))
par.mat$NonDecision <- pars[,"NonDecision"]
par.mat$Separation <- rep(c(pars[,"Separation.2v1"],
pars[,"Separation.3v1"]))
par.mat$Bias <- rep(c(pars[,"Bias.2v1"],
pars[,"Bias.3v1"]))
par.mat$Drift <- c(pars[,"Drift.Control.2v1"],
pars[,"Drift.Control.3v1"],
pars[,"Drift.Inference.2v1"],
pars[,"Drift.Inference.3v1"])
return(par.mat)
}datlist.2v1 <- list(control.2v1.data, inference.2v1.data)
datlist.3v1 <- list(control.3v1.data, inference.3v1.data)
datmat <- list(datlist.2v1,datlist.3v1)
# x has (non-decision, sep-1, bias-1, drift-11, drift-12,
# sep-2, bias-2, drift-21, drift-22)
optim.adult <- optim(p=c(.1, 1,.5, 1, -1, 1, .5, 1, -1), multi.loss.func,
dat=datmat,method="BFGS")
adult.pars <- make.par.mat(optim.adult)
print(adult.pars)## Condition Trial.Type NonDecision Separation Bias Drift
## 1 2v1 Control 0.4170365 1.934698 0.5062048 0.9941540
## 2 3v1 Control 0.4170365 2.501552 0.4413613 1.7742277
## 3 2v1 Inference 0.4170365 1.934698 0.5062048 -0.6624868
## 4 3v1 Inference 0.4170365 2.501552 0.4413613 -1.1155184kid.data <- read.csv('simpimpSCresults_141209_1.csv') %>%
filter(trialtype %in% c("control_double", "inference")) %>%
group_by(subid,trialnum) %>%
mutate(side = if((trialtype == "control_double") & (correct == "Y")) "more"
else if((trialtype == "control_double") & (correct == "N")) "less"
else if((trialtype == "inference") & (correct == "N")) "more"
else "less") %>%
mutate(resp = as.character(factor(side,levels=c("more","less"),
labels=c("upper","lower")))) %>%
mutate(q = reaction.rt/1000) %>%
group_by() %>%
select(q,trialtype,itemnum,resp) %>%
# filter rt, remove rt outside of stdev*2
filter(log(q) < mean(log(q)) + 2 * sd(log(q)),
log(q) > mean(log(q)) - 2 * sd(log(q)))
kid.control.2v1.data <- filter(kid.data,trialtype == "control_double",
itemnum=="2vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(kid.control.2v1.data)
kid.control.3v1.data <- filter(kid.data,trialtype == "control_double",
itemnum=="3vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(kid.control.3v1.data)
kid.inference.2v1.data <- filter(kid.data,trialtype == "inference",
itemnum=="2vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(kid.inference.2v1.data)
kid.inference.3v1.data <- filter(kid.data,trialtype == "inference",
itemnum=="3vs1") %>%
select(q,resp) %>%
as.data.frame
wiener.plot(kid.inference.3v1.data)
Compute kid parameters for each condition separately
kid.optim.control.2v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=kid.control.2v1.data, method="Nelder-Mead")
kid.optim.control.3v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=kid.control.3v1.data, method="Nelder-Mead")
kid.optim.inference.2v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=kid.inference.2v1.data, method="Nelder-Mead")
kid.optim.inference.3v1 <- optim(c(1, .1, .1, 1), wiener_deviance,
dat=kid.inference.3v1.data, method="Nelder-Mead")
kid.indiv.pars <- as.data.frame(rbind(kid.optim.control.2v1$par,
kid.optim.control.3v1$par,
kid.optim.inference.2v1$par,
kid.optim.inference.3v1$par))
names(kid.indiv.pars) <- c("Separation","Non.Decision","Bias","Drift")
kid.indiv.pars$Condition <- c("Control", "Control", "Inference", "Inference")
kid.indiv.pars$Item.Num<- c(2, 3, 2, 3)
print(indiv.pars)## Separation Non.Decision Bias Drift Condition Item.Num
## 1 1.786754 0.4291568 0.4340421 1.1447030 Control 2
## 2 2.996473 0.2842518 0.3047898 2.2676962 Control 3
## 3 1.974870 0.4363610 0.5747067 -0.8593789 Inference 2
## 4 2.632747 0.3087888 0.5847111 -1.4294429 Inference 3print(kid.indiv.pars)## Separation Non.Decision Bias Drift Condition Item.Num
## 1 6.958984 0.5408185 0.6938943 1.2997900 Control 2
## 2 3.626558 0.8303376 0.5382250 0.8859132 Control 3
## 3 2.866273 0.9996504 0.5510521 -0.9332372 Inference 2
## 4 3.035391 1.0323493 0.3808395 -0.5661739 Inference 3Pin parameters for item.num together
# Function for tying parameters together across conditions
kid.datlist.2v1 <- list(kid.control.2v1.data, kid.inference.2v1.data)
kid.datlist.3v1 <- list(kid.control.3v1.data, kid.inference.3v1.data)
# use nlm to estimate parameters
kid.optim.2v1 <- optim(p=c(1, .1, .1, 1, 1), many_drifts,
dat=kid.datlist.2v1,method="Nelder-Mead")
kid.optim.3v1 <- optim(p=c(1, .1, .1, 1, 1), many_drifts,
dat=kid.datlist.3v1,method="Nelder-Mead")
kid.joint.pars <- as.data.frame(rbind(kid.optim.2v1$par,kid.optim.3v1$par))
names(kid.joint.pars) <- c("Separation","Non.Decision","Bias",
"Control", "Inference")
kid.joint.pars <- gather(kid.joint.pars,Condition,Drift,Control:Inference)
kid.joint.pars$Item.Num <- c(2,3,2,3)
kid.joint.pars %<>% select(Separation,Non.Decision,Bias,Drift,
Condition,Item.Num) %>%
arrange(Condition,Item.Num)
print(kid.joint.pars)## Separation Non.Decision Bias Drift Condition Item.Num
## 1 3.685048 0.7149415 0.5837259 1.0472546 Control 2
## 2 3.264180 0.8560350 0.4861925 0.9090952 Control 3
## 3 3.685048 0.7149415 0.5837259 -1.0695535 Inference 2
## 4 3.264180 0.8560350 0.4861925 -0.7468790 Inference 3kid.datlist.2v1 <- list(kid.control.2v1.data, kid.inference.2v1.data)
kid.datlist.3v1 <- list(kid.control.3v1.data, kid.inference.3v1.data)
kid.datmat <- list(kid.datlist.2v1,kid.datlist.3v1)
# x has (non-decision, sep-1, bias-1, drift-11, drift-12,
# sep-2, bias-2, drift-21, drift-22)
optim.kid <- optim(p=c(.1, 1,.5, 1, -1, 1, .5, 1, -1), multi.loss.func,
dat=kid.datmat,method="BFGS")
kid.pars <- make.par.mat(optim.kid)
print(kid.pars)## Condition Trial.Type NonDecision Separation Bias Drift
## 1 2v1 Control 0.7510677 3.564963 0.5934084 1.0142781
## 2 3v1 Control 0.7510677 3.511022 0.4777457 0.9374237
## 3 2v1 Inference 0.7510677 3.564963 0.5934084 -1.0780338
## 4 3v1 Inference 0.7510677 3.511022 0.4777457 -0.7408163print(adult.pars)## Condition Trial.Type NonDecision Separation Bias Drift
## 1 2v1 Control 0.4170365 1.934698 0.5062048 0.9941540
## 2 3v1 Control 0.4170365 2.501552 0.4413613 1.7742277
## 3 2v1 Inference 0.4170365 1.934698 0.5062048 -0.6624868
## 4 3v1 Inference 0.4170365 2.501552 0.4413613 -1.1155184Some copy/pasted code from the paper
#
# set.seed(0)
# dat <- rwiener(n=100, alpha=2, tau=.3, beta=.5, delta=.5)
#
# dat2 <- rwiener(n=100, alpha=2, tau=.3, beta=.5, delta=.75)
#
# dwiener(dat$q[1], alpha=2, tau=.3, beta=.5, delta=.5, resp=dat$resp[1], give_log=FALSE)
#
# curve(dwiener(x, 2, .3, .5, .5, rep("upper", length(x))),
# xlim=c(0,3), main="Density of upper responses",
# ylab="density", xlab="quantile")
#
# pwiener(dat$q[1], alpha=2, tau=.3, beta=.5, delta=.5, resp=dat$resp[1])
#
# wiener_plot(control.data)
#
# x <- c(2, .3, .5, .5)
# wiener_likelihood(x=x, dat=dat)
# wiener_deviance(x=x, dat=control.data)
# wiener_aic(x=x, dat=dat)
# wiener_bic(x=x, dat=dat)
#
# # using optim, first with Nelder-Mead algorithm, then with BFGS
# optim1 <- optim(c(1, .1, .1, 1), wiener_deviance, dat=dat, method="Nelder-Mead")
# optim2 <- optim(optim1[["par"]], wiener_deviance, dat=dat, method="BFGS", hessian=TRUE)
# # using nlm, which uses a Newton-type algorithm
# # create a second data set and a list containing both data sets
# dat2 <- rwiener(n=100, alpha=2, tau=.3, beta=.5, delta=1)
#