chooseOpi("SimHenson")
if(!is.null(opiInitialize(type="C", cap=6)$err))
stop("opiInitialize failed")
#########################################################
# This section is for single location QUESTP
# This example fits a FoS curve
# Note: only fitting threshold and slope,
# modify the domain for FPR and FNR to fit those as well
#########################################################
# 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)
}
#True parameters (variability is determined according to Henson et al. based on threshold)
loc <- list(threshold = 20, fpr = 0.05, fnr = 0.05)
#Function to fit (Gaussian psychometric function)
pSeen <- function(x, params){return(params[3] +
(1 - params[3] - params[4]) *
(1 - pnorm(x, params[1], params[2])))}
#QUEST+
QP <- QUESTP(Fun = pSeen,
stimDomain = list(0:50),
paramDomain = list(seq(0, 40, 1), #Domain for the 50% threshold (Mean)
seq(.5, 8, .5), #Domain for the slope (SD)
seq(0.05, 0.05, 0.05), #Domain for the FPR (static)
seq(0.05, 0.05, 0.05)), #Domain for the FNR (static)
stopType="H", stopValue=4, maxPresentations=500,
makeStim = makeStim,
tt=loc$threshold, fpr=loc$fpr, fnr=loc$fnr,
verbose = 2)
#Plots results
#Henson's FoS function (as implemented in OPI - ground truth)
HensFunction <- function(Th){
SD <- exp(-0.081*Th + 3.27)
SD[SD > 6] <- 6
return(SD)
}
#Stimulus domain
dB_Domain <- 0:50
FoS <- pSeen(dB_Domain, params = QP$final) # Estimated FoS
FoS_GT <- pSeen(dB_Domain, params = c(loc$threshold, HensFunction(loc$threshold),
loc$fpr, loc$fnr)) #Ground truth FoS (based on Henson et al.)
#Plot (seen stimuli at the top, unseen stimuli at the bottom)
plot(dB_Domain, FoS_GT, type = "l", ylim = c(0, 1), xlab = "dB", ylab = "% seen", col = "blue")
lines(dB_Domain, FoS, col = "red")
points(QP$respSeq$stimuli, QP$respSeq$responses, pch = 16,
col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1))
legend("top", inset = c(0, -.2),legend = c("True","Estimated","Stimuli"),
col=c("blue", "red","red"), lty=c(1,1,0),
pch = c(16, 16, 16), pt.cex = c(0, 0, 1),
horiz = TRUE, xpd = TRUE, xjust = 0)
if (!is.null(opiClose()$err))
warning("opiClose() failed")
chooseOpi("SimHenson")
if(!is.null(opiInitialize(type="C", cap=6)$err))
stop("opiInitialize failed")
######################################################################
# This section is for single location QUESTP
# This example shows that QUEST+ can replicate a ZEST procedure
# by fitting a FoS curve with fixed Slope, FPR and FNR
# Compared with ZEST
# Note that QUEST+ should be marginally more efficient in selecting
# the most informative stimulus
######################################################################
# 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)
}
#True parameters (variability is determined according to Henson et al. based on threshold)
loc <- list(threshold = 30, fpr = 0.03, fnr = 0.03)
#Function to fit (Gaussian psychometric function - Fixed slope (same as default likelihood in ZEST))
pSeen <- function(domain, tt){{0.03+(1-0.03-0.03)*(1-pnorm(domain,tt,1))}}
# ZEST-like QUEST+ procedure
QP <- QUESTP(Fun = pSeen,
stimDomain = list(0:40),
paramDomain = list(seq(0, 40, 1)),
stopType="S", stopValue=1.5, maxPresentations=500,
makeStim = makeStim,
tt=loc$threshold, fpr=loc$fpr, fnr=loc$fnr,
verbose = 2)
# ZEST
ZE <- ZEST(domain = 0:40,
stopType="S", stopValue=1.5, maxPresentations=500,
makeStim = makeStim,
tt=loc$threshold, fpr=loc$fpr, fnr=loc$fnr,
verbose = 2)
#Plots results
#Henson's FoS function (as implemented in OPI - ground truth)
HensFunction <- function(Th){
SD <- exp(-0.081*Th + 3.27)
SD[SD > 6] <- 6
return(SD)
}
#Stimulus domain
dB_Domain <- 0:50
FoS_QP <- pSeen(domain = dB_Domain, tt = QP$final) # Estimated FoS
FoS_ZE <- pSeen(domain = dB_Domain, tt = ZE$final) # Estimated FoS
#Plot (seen stimuli at the top, unseen stimuli at the bottom)
plot(dB_Domain, FoS_QP, type = "l", ylim = c(0, 1), xlab = "dB", ylab = "% seen", col = "blue")
lines(dB_Domain, FoS_ZE, col = "red")
points(QP$respSeq$stimuli, QP$respSeq$responses, pch = 16,
col = rgb(red = 0, green = 0, blue = 1, alpha = 0.5))
points(ZE$respSeq[1,], ZE$respSeq[2,], pch = 16,
col = rgb(red = 1, green = 0, blue = 0, alpha = 0.5))
legend("bottomleft", legend = c("QUEST+","ZEST","Stimuli QUEST+", "Stimuli ZEST"),
col=c("blue", "red","blue","red"), lty=c(1,1,0,0),
pch = c(16, 16, 16, 16), pt.cex = c(0, 0, 1, 1),
horiz = FALSE, xpd = TRUE, xjust = 0)
abline(v = loc$threshold, lty = "dashed")
if (!is.null(opiClose()$err))
warning("opiClose() failed")
chooseOpi("SimHenson")
if(!is.null(opiInitialize(type="C", cap=6)$err))
stop("opiInitialize failed")
#########################################################
# This section is for single location QUESTP
# This example fits a broken stick spatial summation function
# with a multi-dimensional stimulus (varying in size and intensity).
# Stimulus sizes are limited to GI, GII, GIII, GIV and GV.
# The example also shows how to use a helper function to
# simulate responses to multi-dimensional stimuli
# (here, the simulated threshold varies based on stimulus size)
#########################################################
makeStim <- function(stim, n) {
s <- list(x=9, y=9, level=dbTocd(stim[1]), size=stim[2], color="white",
duration=200, responseWindow=1500, checkFixationOK=NULL)
class(s) <- "opiStaticStimulus"
return(s)
}
# Helper function for true threshold (depends on log10(stimulus size),
# diameter assumed to be the second element of stim vector)
ttHelper_SS <- function(location) { # returns a function of (stim)
ff <- function(stim) stim
body(ff) <- substitute(
{return(SensF(log10(pi*(stim[2]/2)^2), c(location$Int1, location$Int2, location$Slo2)))}
)
return(ff)
}
# Function of sensivity vs SSize (log10(stimulus area))
SensF <- function(SSize, params){
Sens <- numeric(length(SSize))
for (i in 1:length(SSize)){
Sens[i] <- min(c(params[1] + 10*SSize[i], params[2] + params[3]*SSize[i]))
}
Sens[Sens < 0] <- 0
return(Sens)
}
Sizes <- c(0.1, 0.21, 0.43, 0.86, 1.72)
#True parameters (variability is determined according to Henson et al. based on threshold)
loc <- list(Int1 = 32, Int2 = 28, Slo2 = 2.5, fpr = 0.05, fnr = 0.05, x = 9, y = 9)
# Function to fit (probability of seen given a certain stimulus intensity and size,
# for different parameters)
pSeen <- function(stim, params){
Th <- SensF(log10(pi*(stim[2]/2)^2), params)
return(0.03 +
(1 - 0.03 - 0.03) *
(1 - pnorm(stim[1], Th, 1)))
}
if (FALSE) {
set.seed(111)
#QUEST+ - takes some time to calculate likelihoods
QP <- QUESTP(Fun = pSeen,
stimDomain = list(0:50, Sizes),
paramDomain = list(seq(0, 40, 1), # Domain for total summation intercept
seq(0, 40, 1), # Domain for partial summation intercept
seq(0, 3, 1)), # Domain for partial summation slope
stopType="H", stopValue=1, maxPresentations=500,
makeStim = makeStim,
ttHelper=ttHelper_SS(loc), tt = 30,
fpr=loc$fpr, fnr=loc$fnr,
verbose = 2)
#Stimulus sizes
G <- log10(c(pi*(0.1/2)^2, pi*(0.21/2)^2, pi*(0.43/2)^2, pi*(0.86/2)^2, pi*(1.72/2)^2));
SizesP <- seq(min(G), max(G), .05)
# True and estimated response
Estim_Summation <- SensF(SizesP, params = QP$final) # Estimated spatial summation
GT_Summation <- SensF(SizesP, params = c(loc$Int1, loc$Int2, loc$Slo2)) # True spatial summation
#Plot
plot(10^SizesP, GT_Summation, type = "l", ylim = c(0, 40), log = "x",
xlab = "Stimulus area (deg^2)", ylab = "Sensitivity (dB)", col = "blue")
lines(10^SizesP, Estim_Summation, col = "red")
points(pi*(QP$respSeq$stimuli.2/2)^2, QP$respSeq$stimuli.1, pch = 16,
col = rgb(red = 1, green = 0, blue = 0, alpha = 0.3))
legend("top", inset = c(0, -.2),legend = c("True","Estimated","Stimuli"),
col=c("blue", "red","red"), lty=c(1,1,0),
pch = c(16, 16, 16), pt.cex = c(0, 0, 1),
horiz = TRUE, xpd = TRUE, xjust = 0)
}
if (!is.null(opiClose()$err))
warning("opiClose() failed")
Run the code above in your browser using DataLab