An implementation of the Bayesian test procedures of King-Smith et al.
and Watson and Pelli. Note that we use the term pdf
throughout as in the
original paper, even though they are discrete probability functions in this
implementation.
ZEST(
domain = 0:40,
prior = rep(1/length(domain), length(domain)),
likelihood = sapply(domain, function(tt) 0.03 + (1 - 0.03 - 0.03) * (1 -
stats::pnorm(domain, tt, 1))),
stopType = "S",
stopValue = 1.5,
minStimulus = utils::head(domain, 1),
maxStimulus = utils::tail(domain, 1),
maxSeenLimit = 2,
minNotSeenLimit = 2,
maxPresentations = 100,
minInterStimInterval = NA,
maxInterStimInterval = NA,
verbose = 0,
makeStim,
stimChoice = "mean",
...
)ZEST.start(
domain = 0:40,
prior = rep(1/length(domain), length(domain)),
likelihood = sapply(domain, function(tt) 0.03 + (1 - 0.03 - 0.03) * (1 -
stats::pnorm(domain, tt, 1))),
stopType = "S",
stopValue = 1.5,
minStimulus = utils::head(domain, 1),
maxStimulus = utils::tail(domain, 1),
maxSeenLimit = 2,
minNotSeenLimit = 2,
maxPresentations = 100,
makeStim,
stimChoice = "mean",
...
)
ZEST.step(state, nextStim = NULL, fixedStimValue = NA, fixedResponse = NA)
ZEST.stop(state)
ZEST.final(state)
ZEST
returns a list containing
npres
Total number of presentations used.
respSeq
Response sequence stored as a matrix: row 1 is dB values of stimuli, row 2
is 1/0 for seen/not-seen, row 3 is fixated 1/0 (always 1 if checkFixationOK
not
present in stim objects returned from makeStim
).
pdfs
If verbose
is bigger than 0, then this is a list of the pdfs used for each presentation, otherwise NULL.
final
The mean/median/mode of the final pdf, depending on stimChoice
, which is the determined threshold.
opiResp
A list of responses received from each successful call to opiPresent
within ZEST
.
ZEST.start
returns a list that can be passed to ZEST.step
, ZEST.stop
, and
ZEST.final
. It represents the state of a ZEST at a single location at a point in time
and contains the following.
name
ZEST
.
pdf
Current pdf: vector of probabilities the same length as domain
.
numPresentations
The number of times ZEST.step
has been called on this state.
stimuli
A vector containing the stimuli used at each call of ZEST.step
.
responses
A vector containing the responses received at each call of ZEST.step
.
responseTimes
A vector containing the response times received at each call of ZEST.step
.
fixated
A vector containing TRUE/FALSE if fixation was OK according to
checkFixationOK
for each call of ZEST.step
(defaults to TRUE if
checkFixationOK
not present).
opiResp
A list of responses received from each call to opiPresent
within ZEST.step
.
A copy of all of the parameters supplied to ZEST.start: domain
,
likelihood
, stopType
, stopValue
, minStimulus
, maxStimulus
,
maxSeenLimit
, minNotSeenLimit
, maxPresentations
, makeStim
,
stimChoice
, currSeenLimit
, currNotSeenLimit
, and opiParams
.
ZEST.step
returns a list containing
* state
The new state after presenting a stimuli and getting a response.
* resp
The return from the opiPresent
call that was made.
ZEST.stop
returns TRUE
if the ZEST has reached its stopping criteria, and FALSE
otherwise.
ZEST.final
returns an estimate of threshold based on state. If state$stimChoice
is mean
then the mean is returned. If state$stimChoice
is mode
then the
mode is returned. If state$stimChoice
is median
then the median is returned.
A list containing
* state
The new state after presenting a stimuli and getting a response.
* resp
The return from the opiPresent
call that was made.
Vector of values over which pdf is kept.
Starting probability distribution over domain. Same length as domain
.
Matrix where likelihood[s,t]
is likelihood of seeing s
given t
is the true threshold. That is, Pr(s|t) where s
and t
are
indexes into domain
.
N
, for number of presentations; S
, for standard deviation
of the pdf; and H
, for the entropy of the pdf.
Value for number of presentations (stopType=N
), standard deviation
(stopType=S)
or Entropy (stopType=H
).
The smallest stimuli that will be presented. Could be different from
domain[1]
.
The largest stimuli that will be presented. Could be different from
tail(domain,1)
.
Will terminate if maxStimulus
value is seen this many times.
Will terminate if minStimulus
value is not seen this many times.
Maximum number of presentations regardless of stopType
.
If both minInterStimInterval
and maxInterStimInterval
are not NA
, then between each stimuli there is a random wait period drawn uniformly
between minInterStimInterval
and maxInterStimInterval
.
minInterStimInterval
.
verbose=0
does nothing, verbose=1
stores pdfs for returning,
and verbose=2
stores pdfs and also prints each presentation.
A function that takes a dB value and numPresentations and returns an OPI datatype ready for passing to opiPresent. See examples.
A true ZEST procedure uses the "mean"
of the current pdf as the stimulus,
but "median"
and "mode"
(as used in a QUEST procedure) are provided for your
enjoyment.
Extra parameters to pass to the opiPresent function
Current state of the ZEST returned by ZEST.start
and ZEST.step
.
A valid object for opiPresent
to use as its nextStim
.
A number in state$domain
that, is !is.na
, will be used as the stimulus value
overriding state$minStimulus
, state$maxStimulus
and state$stimChoice
.
Ignored if !is.na
otherwise used as the response to the stim shown.
This is an implementation of King-Smith et al.'s ZEST procedure and Watson and Pelli's QUEST procedure. All presentations are rounded to an element of the supplied domain.
Note this function will repeatedly call opiPresent
for a stimulus until opiPresent
returns NULL
(ie no error occurred).
The checkFixationOK
function is called (if present in stim made from makeStim
)
after each presentation, and if it returns FALSE, the pdf for that location is not changed
(ie the presentation is ignored), but the stim, number of presentations etc is recorded in
the state.
If more than one ZEST is to be interleaved (for example, testing multiple locations), then the
ZEST.start
, ZEST.step
, ZEST.stop
and ZEST.final
calls can maintain
the state of the ZEST after each presentation, and should be used. If only a single ZEST is
required, then the simpler ZEST
can be used, which is a wrapper for the four functions
that maintain state. See examples below.
P.E. King-Smith, S.S. Grigsny, A.J. Vingrys, S.C. Benes, and A. Supowit. "Efficient and Unbiased Modifications of the QUEST Threshold Method: Theory, Simulations, Experimental Evaluation and Practical Implementation", Vision Research 34(7) 1994. Pages 885-912.
A.B. Watson and D.G. Pelli. "QUEST: A Bayesian adaptive psychophysical method", Perception and Psychophysics 33 1983. Pages 113-l20.
A. Turpin, P.H. Artes and A.M. McKendrick "The Open Perimetry Interface: An enabling tool for clinical visual psychophysics", Journal of Vision 12(11) 2012.
dbTocd
, opiPresent
chooseOpi("SimHenson")
if(!is.null(opiInitialize(type="C", cap=6)$err))
stop("opiInitialize failed")
##############################################
# This section is for single location ZESTs
##############################################
# Stimulus is Size III white-on-white as in the HFA
makeStim <- function(db, n) {
s <- list(x=9, y=9, level=dbTocd(db), size=0.43, color="white",
duration=200, responseWindow=1500, checkFixationOK=NULL)
class(s) <- "opiStaticStimulus"
return(s)
}
repp <- function(...) sapply(1:50, function(i) ZEST(makeStim=makeStim, ...))
a <- repp(stopType="H", stopValue= 3, verbose=0, tt=30, fpr=0.03)
b <- repp(stopType="S", stopValue=1.5, verbose=0, tt=30, fpr=0.03)
c <- repp(stopType="S", stopValue=2.0, verbose=0, tt=30, fpr=0.03)
d <- repp(stopType="N", stopValue= 50, verbose=0, tt=30, fpr=0.03)
e <- repp(prior=dnorm(0:40,m=0,s=5), tt=30, fpr=0.03)
f <- repp(prior=dnorm(0:40,m=10,s=5), tt=30, fpr=0.03)
g <- repp(prior=dnorm(0:40,m=20,s=5), tt=30, fpr=0.03)
h <- repp(prior=dnorm(0:40,m=30,s=5), tt=30, fpr=0.03)
layout(matrix(1:2,1,2))
boxplot(lapply(list(a,b,c,d,e,f,g,h), function(x) unlist(x["final",])))
boxplot(lapply(list(a,b,c,d,e,f,g,h), function(x) unlist(x["npres",])))
##############################################
# This section is for multiple ZESTs
##############################################
makeStimHelper <- function(db,n, x, y) { # returns a function of (db,n)
ff <- function(db, n) db+n
body(ff) <- substitute({
s <- list(x=x, y=y, level=dbTocd(db), size=0.43, color="white",
duration=200, responseWindow=1500, checkFixationOK=NULL)
class(s) <- "opiStaticStimulus"
return(s)
}, list(x=x,y=y))
return(ff)
}
# List of (x, y, true threshold) triples
locations <- list(c(9,9,30), c(-9,-9,32), c(9,-9,31), c(-9,9,33))
# Setup starting states for each location
states <- lapply(locations, function(loc) {
ZEST.start(
domain=-5:45,
minStimulus=0,
maxStimulus=40,
makeStim=makeStimHelper(db,n,loc[1],loc[2]),
stopType="S", stopValue= 1.5, tt=loc[3], fpr=0.03, fnr=0.01)})
# Loop through until all states are "stop"
while(!all(st <- unlist(lapply(states, ZEST.stop)))) {
i <- which(!st) # choose a random,
i <- i[runif(1, min=1, max=length(i))] # unstopped state
r <- ZEST.step(states[[i]]) # step it
states[[i]] <- r$state # update the states
}
finals <- lapply(states, ZEST.final) # get final estimates of threshold
for(i in 1:length(locations)) {
#cat(sprintf("Location (%+2d,%+2d) ",locations[[i]][1], locations[[i]][2]))
#cat(sprintf("has threshold %4.2f\n", finals[[i]]))
}
if (!is.null(opiClose()$err))
warning("opiClose() failed")
Run the code above in your browser using DataLab