Title: | Functions for Analyzing Psychophysical Data in R |
---|---|
Description: | An assortment of functions that could be useful in analyzing data from psychophysical experiments. It includes functions for calculating d' from several different experimental designs, links for m-alternative forced-choice (mafc) data to be used with the binomial family in glm (and possibly other contexts) and self-Start functions for estimating gamma values for CRT screen calibrations. |
Authors: | Kenneth Knoblauch |
Maintainer: | Ken Knoblauch <[email protected]> |
License: | GPL |
Version: | 0.3 |
Built: | 2024-12-11 07:01:28 UTC |
Source: | CRAN |
An assortment of functions that could be useful in analyzing data from pyschophysical experiments. It includes functions for calculating d' from several different experimental designs, links for mafc to be used with the binomial family in glm (and possibly other contexts) and a self-Start function for estimating gamma values for CRT screen calibrations.
For the moment, the package contains several functions for calculating
for a variety of psychophysical paradigms, some link functions for
the binomial family in glm (and perhaps other functions) for fitting
psychometric functions from mAFC experiments and a self-Start function
for estimating the value of the exponent, gamma, based on the
luminance calibration of the three channels of a CRT-like display.
Kenneth Knoblauch <[email protected]>
Calulate for ABX paradigm either
assuming a differencing strategy or independent observations
dprime.ABX(Hits, FA, zdiff, Pc.unb, method = "diff")
dprime.ABX(Hits, FA, zdiff, Pc.unb, method = "diff")
Hits |
numeric in [0, 1] corresponding to Hit rate |
FA |
numeric in [0, 1] corresponding to False alarm rate |
zdiff |
numeric. Difference of z-scores for Hit and False Alarm rates |
Pc.unb |
numeric in [0, 1]. Proportion correct for an unbiased observer,
|
method |
character. Specifies the model to describe the observer's criterion for dividing up the decision space, must be either "diff" for a differencing strategy (the default) or "IO" for independent observations. |
Two different strategies have been described for how the
observer partitions the decision space in the ABX
paradigm, either based on Independent Observations of
each stimulus or on a differencing strategy. The differecing
strategy is the default. can be calculated
either from the
H
and FA
rates, from the difference of
z-scores or from the probability correct of an unbiased observer.
Returns the value of
Kenneth Knoblauch
MacMillan, N. A. and Creeman, C. D. (1991) Detection Theory: A User's Guide Cambridge University Press
Green, D. M. and Swets, J. A. (1966) Signal Detection Theory and Psychophysics Robert E. Krieger Publishing Company
dprime.mAFC
, dprime.SD
,
dprime.oddity
dprime.ABX(H = 0.75, F = 0.3) dprime.ABX(H = 0.75, F = 0.3, method = "IO") dprime.ABX(zdiff = qnorm(0.75) - qnorm(0.3)) dprime.ABX(Pc.unb = pnorm( (qnorm(0.75) - qnorm(0.3))/2 ))
dprime.ABX(H = 0.75, F = 0.3) dprime.ABX(H = 0.75, F = 0.3, method = "IO") dprime.ABX(zdiff = qnorm(0.75) - qnorm(0.3)) dprime.ABX(Pc.unb = pnorm( (qnorm(0.75) - qnorm(0.3))/2 ))
Calculate the value of for an m-alternative forced choice paradigm
dprime.mAFC(Pc, m)
dprime.mAFC(Pc, m)
Pc |
The proportion of correct responses based on either the Hit rate or based on an unbiased observer |
m |
The number of alternative choices, an integer > 1. |
The probability of a correct response in m-alternative forced-choice, assuming independence, is based on the product of the likelihoods of the signal alternative generating the strongest response and the m - 1 noise alternatives generating responses less than this value (Green and Dai, 1991). For a Yes-No paradigm, the sensitivity is calculated more simply as
where and
are the Hit and False Alarm rates,
respectively.
Returns the value of
Currently is only valid for in the interval [-10, 10]
which should be well outside the range of sensory differences that
this paradigm is used to investigate.
Kenneth Knoblauch
Green, D. M. and Dai, H. (1991) Probability of being correct with 1 of M orthogonal signals. Perception & Psychophysics, 49, 100–101.
Green, D. M. and Swets, J. A. (1966) Signal Detection Theory and Psychophysics Robert E. Krieger Publishing Company
See Also dprime.ABX
, dprime.SD
,
dprime.oddity
dprime.mAFC(0.8, 4)
dprime.mAFC(0.8, 4)
Calculate for a 3 stimulus (triangular) paradigm. Two of the
stimuli are the same and the observer must designate the stimulus that
is different.
dprime.oddity(Pc.tri)
dprime.oddity(Pc.tri)
Pc.tri |
numeric in (1/3, 1). The proportion of correct responses for an unbiased observer. |
Returns the value of
Kenneth Knoblauch
Frijters, G. S., Kooistra, A. and Verijken, P. F. G. (1980) Tables of
for the triangular method and the 3-AFC signal detection procedure.
Perception & Psychophysics, 27, 176–178.
MacMillan, N. A. and Creeman, C. D. (1991) Detection Theory: A User's Guide Cambridge University Press
Green, D. M. and Swets, J. A. (1966) Signal Detection Theory and Psychophysics Robert E. Krieger Publishing Company
dprime.mAFC
, dprime.SD
,
dprime.ABX
dprime.oddity(0.8)
dprime.oddity(0.8)
Calulate for same-different paradigm either
assuming a differencing strategy or independent observations
dprime.SD(H, FA, zdiff, Pcmax, method = "diff")
dprime.SD(H, FA, zdiff, Pcmax, method = "diff")
H |
numeric in [0, 1] corresponding to Hit rate |
FA |
numeric in [0, 1] corresponding to False alarm rate |
zdiff |
numeric. Difference of z-scores for Hit and False Alarm rates ( only valid for method "IO") |
Pcmax |
numeric in [0, 1]. Proportion correct for an unbiased observer,
|
method |
character. Specifies the model to describe the observer's criterion for dividing up the decision space, must be either "diff" for a differencing strategy (the default) or "IO" for independent observations. |
Two different strategies have been described for how the
observer partitions the decision space in the same-different
paradigm. With Independent Observations, can be calculated
either from the
H
and FA
rates, from the difference of
z-scores or from the probability correct of an unbiased observer.
Only one of these three choices should be specified in the arguments.
For the differencing strategy, only the first of these choices is valid.
Returns the value of
Kenneth Knoblauch
MacMillan, N. A. and Creeman, C. D. (1991) Detection Theory: A User's Guide Cambridge University Press
Green, D. M. and Swets, J. A. (1966) Signal Detection Theory and Psychophysics Robert E. Krieger Publishing Company
dprime.mAFC
, dprime.ABX
,
dprime.oddity
dprime.SD(H = 0.642, F = 0.3) dprime.SD(H = 0.75, F = 0.3, method = "IO") dprime.SD(zdiff = qnorm(0.75) - qnorm(0.3), method = "IO") dprime.SD(Pcmax = pnorm( (qnorm(0.75) - qnorm(0.3))/2 ), method = "IO")
dprime.SD(H = 0.642, F = 0.3) dprime.SD(H = 0.75, F = 0.3, method = "IO") dprime.SD(zdiff = qnorm(0.75) - qnorm(0.3), method = "IO") dprime.SD(Pcmax = pnorm( (qnorm(0.75) - qnorm(0.3))/2 ), method = "IO")
Letter detection and identification at 2 degrees eccentricity in the visual field. On each trial, one of four letters (b, d, p, q) were presented in one of four positions (superior, inferior, left, right) in the visual field. In a given session, the letter height was fixed. Six contrast levels were tested in each session. The data indicate the proportion of correctly identified positions, referred to here as detection, and the proportion of correctly identified letters, conditional on correct identification.
data(ecc2)
data(ecc2)
A data frame with 48 observations on the following 5 variables.
Contr
numeric. The contrast of the stimulus, defined as Weberian contrast.
task
a factor with levels DET
ID
indicating the two tasks, detection and identification.
Size
a numeric vector indicating the letter height
Correct
an integer vector indicating the number of correct responses (DET
or ID
).
Incorrect
an integer vector, indicating the number of incorrect responses.
Yssaad-Fesselier, R. and Knoblauch, K. (2006) Modeling psychometric functions in R. Behav Res Methods., 38(1), 28–41.
data(ecc2) library(lattice) xyplot(Correct/(Correct + Incorrect) ~ Contr | Size * task, ecc2, type = "b", scale = list(x = list(log = TRUE), y = list(limits = c(0, 1.05))), xlab = "Contrast", ylab = "Proportion Correct Response", panel = function(x, y, ...) { panel.xyplot(x, y, ...) panel.abline(h = 0.25, lty = 2) })
data(ecc2) library(lattice) xyplot(Correct/(Correct + Incorrect) ~ Contr | Size * task, ecc2, type = "b", scale = list(x = list(log = TRUE), y = list(limits = c(0, 1.05))), xlab = "Contrast", ylab = "Proportion Correct Response", panel = function(x, y, ...) { panel.xyplot(x, y, ...) panel.abline(h = 0.25, lty = 2) })
A wrapper for glm
in which the deviance for the model with binomial family and link probit.lambda
is profiled as a function of lambda
, the upper asymptote of the psychometric function.
glm.lambda(formula, data, NumAlt = 2, lambda = seq(0, 0.1, len = 40), plot.it = FALSE, ...)
glm.lambda(formula, data, NumAlt = 2, lambda = seq(0, 0.1, len = 40), plot.it = FALSE, ...)
formula |
a symbolic description of the model to be fit |
data |
an optional data frame, list or environment (or object coercible by |
NumAlt |
the number of alternatives, |
lambda |
a sequence of values to profile for the upper asymptote of the psychometric function |
plot.it |
logical indicating whether to plot the profile of the deviances as a function of |
... |
further arguments passed to |
The psychometric function fit to the data is described by
where is the number of alternatives and the lower asymptote,
is the upper asymptote and
is the cumulative normal function.
returns an object of class ‘lambda’ which inherits from classes ‘glm’ and ‘lm’. It only differs from an object of class ‘glm’ in including two additional components, lambda
, giving the estimated minimum of the profile by fitting a quadratic to the profile and a data frame containing the profiled deviance values for each value of lambda
tested. The degrees of freedom are reduced by 1 to take into account the estimation of lambda
.
If the minimum occurs outside the interval examined, an error might occur. In any case, re-running the function with a new range of lambda
that includes the minimum should work. if the plotted profile indicates that the fitted quadratic does not describe well the profile at the minimum, refitting with a more restricted range of lambda
is recommended.
Ken Knoblauch
Wichmann, F. A. and Hill, N. J. (2001) The psychometric function: I.Fitting, sampling, and goodness of fit. Percept Psychophys., 63(8), 1293–1313.
Yssaad-Fesselier, R. and Knoblauch, K. (2006) Modeling psychometric
functions in R. Behav Res Methods., 38(1), 28–41. (for examples
with gnlr
).
mafc
, glm
, probit.lambda
, family
b <- 3.5 g <- 1/3 d <- 0.025 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses set.seed(12161952) truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ## First with upper asymptote at 1 dd.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.probit(3))) summary(dd.glm) dd.lam <- glm.lambda(resp.mat ~ cnt, NumAlt = 3, lambda = seq(0, 0.03, len = 100), plot.it = TRUE) summary(dd.lam) ## can refine interval, but doesn't change result much dd.lam2 <- glm.lambda(resp.mat ~ cnt, NumAlt = 3, lambda = seq(dd.lam$lambda/sqrt(2), dd.lam$lambda*sqrt(2), len = 100), plot.it = TRUE) summary(dd.lam2) ## Compare fits w/ and w/out lambda anova(dd.glm, dd.lam2, test = "Chisq") plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(dd.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) lines(pcnt, predict(dd.lam, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2)
b <- 3.5 g <- 1/3 d <- 0.025 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses set.seed(12161952) truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ## First with upper asymptote at 1 dd.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.probit(3))) summary(dd.glm) dd.lam <- glm.lambda(resp.mat ~ cnt, NumAlt = 3, lambda = seq(0, 0.03, len = 100), plot.it = TRUE) summary(dd.lam) ## can refine interval, but doesn't change result much dd.lam2 <- glm.lambda(resp.mat ~ cnt, NumAlt = 3, lambda = seq(dd.lam$lambda/sqrt(2), dd.lam$lambda*sqrt(2), len = 100), plot.it = TRUE) summary(dd.lam2) ## Compare fits w/ and w/out lambda anova(dd.glm, dd.lam2, test = "Chisq") plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(dd.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) lines(pcnt, predict(dd.lam, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2)
A probit fit of a psychometric function with upper asymptote less than 1 is obtained by cycling between a fit with glm
using the probit.lambda
link and optimize
to estimate lambda
, 1 - the upper asymptotic value, until the log Likelihood changes by less than a pre-set tolerance.
glm.WH(formula, data, NumAlt = 2, lambda.init = 0.01, interval = c(0, 0.05), trace = FALSE, tol = 1e-06, ...)
glm.WH(formula, data, NumAlt = 2, lambda.init = 0.01, interval = c(0, 0.05), trace = FALSE, tol = 1e-06, ...)
formula |
a symbolic description of the model to be fit. |
data |
an optional data frame, list or enviroment (or object coercible by |
NumAlt |
integer indicating the number of alternatives (> 1) in the mafc-task. (Default: 2). |
lambda.init |
numeric, initial estimate of 1 - upper asymptote. |
interval |
numeric vector giving interval endpoints within which to search for |
trace |
logical, indicating whether or not to print out a trace of the iterative process. |
tol |
numeric, tolerance for ending iterations. |
... |
futher arguments passed to |
The psychometric function fit to the data is described by
where is the number of alternatives and the lower asymptote,
is the upper asymptote and
is the cumulative normal function.
returns an object of class ‘lambda’ which inherits from classes ‘glm’ and ‘lm’. It only differs from an object of class ‘glm’ in including an additional components, lambda
, giving the estimated minimum of lambda
. The degrees of freedom are reduced by 1 to take into account the estimation of lambda
.
Ken Knoblauch
Wichmann, F. A. and Hill, N. J. (2001) The psychometric function: I.Fitting, sampling, and goodness of fit. Percept Psychophys., 63(8), 1293–1313.
Yssaad-Fesselier, R. and Knoblauch, K. (2006) Modeling psychometric
functions in R. Behav Res Methods., 38(1), 28–41. (for examples
with gnlr
).
mafc
, glm
,glm.lambda
, probit.lambda
, family
b <- 3.5 g <- 1/4 d <- 0.04 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) tst.glm <- glm(resp.mat ~ cnt, binomial(mafc.probit(1/g))) pcnt <- seq(0.005, 1, len = 1000) plot(cnt, phat, log = "x", ylim = c(0, 1), xlim = c(0.005, 1), cex = 1.75) lines(pcnt, predict(tst.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) tst.lam <- glm.WH(resp.mat ~ cnt, NumAlt = 1/g, trace = TRUE) lines(pcnt, predict(tst.lam, data.frame(cnt = pcnt), type = "response"), lty = 2, lwd = 2)
b <- 3.5 g <- 1/4 d <- 0.04 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) tst.glm <- glm(resp.mat ~ cnt, binomial(mafc.probit(1/g))) pcnt <- seq(0.005, 1, len = 1000) plot(cnt, phat, log = "x", ylim = c(0, 1), xlim = c(0.005, 1), cex = 1.75) lines(pcnt, predict(tst.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) tst.lam <- glm.WH(resp.mat ~ cnt, NumAlt = 1/g, trace = TRUE) lines(pcnt, predict(tst.lam, data.frame(cnt = pcnt), type = "response"), lty = 2, lwd = 2)
These functions provide links for the binamial family so that psychometric functions can be fit with both the upper and lower asymptotes different from 1 and 0, respectively.
logit.2asym(g, lam) probit.2asym(g, lam) cauchit.2asym(g, lam) cloglog.2asym(g, lam) weib.2asym( ... )
logit.2asym(g, lam) probit.2asym(g, lam) cauchit.2asym(g, lam) cloglog.2asym(g, lam) weib.2asym( ... )
g |
numeric in the range (0, 1), normally <= 0.5, however, which specifies the lower asymptote of the psychometric function. |
lam |
numeric in the range (0, 1), specifying 1 - the upper asymptote of the psychometric function. |
... |
used just to pass along the formals of |
These links are used to specify psychometric functions with the form
where is the lower asymptote and
is
the upper asymptote, and
is the base psychometric function, varying between 0 and 1.
Each link returns a list containing functions required for relating the response to the linear predictor in generalized linear models and the name of the link.
linkfun |
The link function |
linkinv |
The inverse link function |
mu.eta |
The derivative of the inverse link |
valideta |
The domain over which the linear predictor is valid |
link |
A name to be used for the link |
Kenneth Knoblauch
Klein S. A. (2001) Measuring, estimating, and understanding the psychometric function: a commentary. Percept Psychophys., 63(8), 1421–1455.
Wichmann, F. A. and Hill, N. J. (2001) The psychometric function: I.Fitting, sampling, and goodness of fit. Percept Psychophys., 63(8), 1293–1313.
glm
, glm
make.link
, psyfun.2asym
#A toy example, b <- 3 g <- 0.05 # simulated false alarm rate d <- 0.03 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- psyfun.2asym(resp.mat ~ cnt, link = probit.2asym) ddlog.glm <- psyfun.2asym(resp.mat ~ cnt, link = logit.2asym) # Can fit a Weibull function, but use log contrast as variable ddweib.glm <- psyfun.2asym(resp.mat ~ log(cnt), link = weib.2asym) ddcau.glm <- psyfun.2asym(resp.mat ~ cnt, link = cauchit.2asym) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 5) lines(pcnt, predict(ddlog.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2, col = "blue") lines(pcnt, predict(ddweib.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey") lines(pcnt, predict(ddcau.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey", lty = 2)
#A toy example, b <- 3 g <- 0.05 # simulated false alarm rate d <- 0.03 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- psyfun.2asym(resp.mat ~ cnt, link = probit.2asym) ddlog.glm <- psyfun.2asym(resp.mat ~ cnt, link = logit.2asym) # Can fit a Weibull function, but use log contrast as variable ddweib.glm <- psyfun.2asym(resp.mat ~ log(cnt), link = weib.2asym) ddcau.glm <- psyfun.2asym(resp.mat ~ cnt, link = cauchit.2asym) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 5) lines(pcnt, predict(ddlog.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2, col = "blue") lines(pcnt, predict(ddweib.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey") lines(pcnt, predict(ddcau.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey", lty = 2)
These provide links for the binomial family for fitting m-alternative forced-choice psychophysical functions.
mafc.logit( .m = 2 ) mafc.probit( .m = 2 ) mafc.cloglog( .m = 2 ) mafc.weib( ... ) mafc.cauchit( .m = 2 )
mafc.logit( .m = 2 ) mafc.probit( .m = 2 ) mafc.cloglog( .m = 2 ) mafc.weib( ... ) mafc.cauchit( .m = 2 )
.m |
is the integer number (>1) of choices (Default to 2AFC). For m = 1 (Yes/No paradigm), use one of the built-in links for the binomial family. |
... |
just to pass along the formals of |
These
functions provide links for fitting psychometric functions arising from
an m-alternative forced-choice experiment. The estimated coefficients
of the linear predictor influence both the location and the slope of the
psychometric function(s), but provide no means of estimating the upper
aymptote which is constrained to approach 1. If the upper asympotote
must be estimated, it would be better to maximize directly the
likelihood, either with a function like optim
or gnlr
from
package gnlm (available at
https://www.commanster.eu/rcode.html). Alternatively,
the function probit.lambda
can be used with a known
upper asymptote, or glm.lambda
or glm.WH
to estimate one, with a probit link. mafc.weib
is just an
alias for mafc.cloglog
.
Each link returns a list containing functions required for relating the response to the linear predictor in generalized linear models and the name of the link.
linkfun |
The link function |
linkinv |
The inverse link function |
mu.eta |
The derivative of the inverse link |
valideta |
The domain over which the linear predictor is valid |
link |
A name to be used for the link |
Kenneth Knoblauch
Williams J, Ramaswamy D and Oulhaj A (2006) 10 Hz flicker improves recognition memory in older people BMC Neurosci. 2006 5;7:21 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1434755/ (for an example developed prior to this one, but for m = 2).
Klein S. A. (2001) Measuring, estimating, and understanding the psychometric function: a commentary. Percept Psychophys., 63(8), 1421–1455.
Wichmann, F. A. and Hill, N. J. (2001) The psychometric function: I.Fitting, sampling, and goodness of fit. Percept Psychophys., 63(8), 1293–1313.
Yssaad-Fesselier, R. and Knoblauch, K. (2006) Modeling psychometric
functions in R. Behav Res Methods., 38(1), 28–41. (for examples
with gnlr
).
family
, make.link
, glm
,
optim
, probit.lambda
, glm.lambda
,
glm.WH
#A toy example, b <- 3.5 g <- 1/3 d <- 0.0 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.probit(3))) ddlog.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.logit(3))) # Can fit a Weibull function, but use log contrast as variable ddweib.glm <- glm(resp.mat ~ log(cnt), family = binomial(mafc.cloglog(3))) ddcau.glm <- glm(resp.mat ~ log(cnt), family = binomial(mafc.cauchit(3))) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) lines(pcnt, predict(ddlog.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2) lines(pcnt, predict(ddweib.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey") lines(pcnt, predict(ddcau.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey", lty = 2) # Weibull parameters \alpha and \beta cc <- coef(ddweib.glm) alph <- exp(-cc[1]/cc[2]) bet <- cc[2] #More interesting example with data from Yssaad-Fesselier and Knoblauch data(ecc2) ecc2.glm <- glm(cbind(Correct, Incorrect) ~ Contr * Size * task, family = binomial(mafc.probit(4)), data = ecc2) summary(ecc2.glm) ecc2$fit <- fitted(ecc2.glm) library(lattice) xyplot(Correct/(Correct + Incorrect) ~ Contr | Size * task, data = ecc2, subscripts = TRUE, ID = with(ecc2, Size + as.numeric(task)), scale = list(x = list(log = TRUE), y = list(limits = c(0, 1.05))), xlab = "Contrast", ylab = "Proportion Correct Response", aspect = "xy", panel = function(x, y, subscripts, ID, ...) { which = unique(ID[subscripts]) llines(x, ecc2$fit[which ==ID], col = "black", ...) panel.xyplot(x, y, pch = 16, ...) panel.abline(h = 0.25, lty = 2, ...) } )
#A toy example, b <- 3.5 g <- 1/3 d <- 0.0 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.probit(3))) ddlog.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.logit(3))) # Can fit a Weibull function, but use log contrast as variable ddweib.glm <- glm(resp.mat ~ log(cnt), family = binomial(mafc.cloglog(3))) ddcau.glm <- glm(resp.mat ~ log(cnt), family = binomial(mafc.cauchit(3))) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) lines(pcnt, predict(ddlog.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2) lines(pcnt, predict(ddweib.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey") lines(pcnt, predict(ddcau.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey", lty = 2) # Weibull parameters \alpha and \beta cc <- coef(ddweib.glm) alph <- exp(-cc[1]/cc[2]) bet <- cc[2] #More interesting example with data from Yssaad-Fesselier and Knoblauch data(ecc2) ecc2.glm <- glm(cbind(Correct, Incorrect) ~ Contr * Size * task, family = binomial(mafc.probit(4)), data = ecc2) summary(ecc2.glm) ecc2$fit <- fitted(ecc2.glm) library(lattice) xyplot(Correct/(Correct + Incorrect) ~ Contr | Size * task, data = ecc2, subscripts = TRUE, ID = with(ecc2, Size + as.numeric(task)), scale = list(x = list(log = TRUE), y = list(limits = c(0, 1.05))), xlab = "Contrast", ylab = "Proportion Correct Response", aspect = "xy", panel = function(x, y, subscripts, ID, ...) { which = unique(ID[subscripts]) llines(x, ecc2$fit[which ==ID], col = "black", ...) panel.xyplot(x, y, pch = 16, ...) panel.abline(h = 0.25, lty = 2, ...) } )
This function provides a link for the binomial family for fitting m-alternative forced-choice, with a probit link and with the upper asymptote permitted to be less than 1.
probit.lambda(m = 2, lambda = 0)
probit.lambda(m = 2, lambda = 0)
m |
is the integer number (>1) of choices (Default to 2AFC). |
lambda |
number in [0, 1] indicating 1 minus the upper asymptotic value of the psychometric function. |
This function provides a link for fitting psychometric functions arising from an m-alternative forced-choice experiment using a probit link and allowing that the upper aymptote is less than 1. The psychometric function fit to the data is described by
where is the number of alternatives and the lower asymptote,
is the upper asymptote and
is the cumulative normal function.
The link returns a list containing functions required for relating the response to the linear predictor in generalized linear models and the name of the link.
linkfun |
The link function |
linkinv |
DTHe inverse link function |
mu.eta |
The derivative of the inverse link function |
valideta |
The domain over which the linear predictor is valid |
link |
A name to be used for the link |
Due to the difficulty of the task, subject error or incorrectly recorded data, psychophysical data may reveal less than perfect performance when stimulus differences are readily visible. When this occurs, letting the upper asymptote be less than 1 often results in a better fit to the data and a less-biased estimate of the steepness of the curve (see example below).
Ken Knoblauch
Wichmann, F. A. and Hill, N. J. (2001) The psychometric function: I.Fitting, sampling, and goodness of fit. Percept Psychophys., 63(8), 1293–1313.
mafc
, glm
, glm.lambda
, family
, make.link
b <- 3.5 g <- 1/3 d <- 0.025 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.probit(3))) ddprob.lam <- glm(resp.mat ~ cnt, family = binomial(probit.lambda(3, 0.025))) AIC(ddprob.glm, ddprob.lam) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) lines(pcnt, predict(ddprob.lam, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2)
b <- 3.5 g <- 1/3 d <- 0.025 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- glm(resp.mat ~ cnt, family = binomial(mafc.probit(3))) ddprob.lam <- glm(resp.mat ~ cnt, family = binomial(probit.lambda(3, 0.025))) AIC(ddprob.glm, ddprob.lam) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2) lines(pcnt, predict(ddprob.lam, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2)
Fits psychometric functions allowing for variation of both upper and lower asymptotes. Uses a procedure that alternates between fitting linear predictor with glm
and estimating the asymptotes with optim
until a minimum in -log likelihood is obtained within a tolerance.
psyfun.2asym(formula, data, link = logit.2asym, init.g = 0.01, init.lam = 0.01, trace = FALSE, tol = 1e-06, mxNumAlt = 50, ...)
psyfun.2asym(formula, data, link = logit.2asym, init.g = 0.01, init.lam = 0.01, trace = FALSE, tol = 1e-06, mxNumAlt = 50, ...)
formula |
a two sided formula specifying the response and the linear predictor |
data |
a data frame within which the formula terms are interpreted |
link |
a link function for the binomial family that allows specifying both upper and lower asymptotes |
init.g |
numeric specifying the initial estimate for the lower asymptote |
init.lam |
numeric specifying initial estimate for 1 - upper asymptote |
trace |
logical indicating whether to show the trace of the minimization of -log likelihood |
tol |
numeric indicating change in -log likelihood as a criterion for stopping iteration. |
mxNumAlt |
integer indicating maximum number of alternations between |
... |
additional arguments passed to |
The function is a wrapper for glm
for fitting psychometric functions with the equation
where is the lower asymptote and
is
the upper asymptote, and
is the base psychometric function, varying between 0 and 1.
list of class ‘lambda’ inheriting from classes ‘glm’ and ‘lm’ and containing additional components
lambda |
numeric indicating 1 - upper asymptote |
gam |
numeric indicating lower asymptote |
SElambda |
numeric indicating standard error estimate for lambda based on the Hessian of the last interation of |
SEgam |
numeric indicating standard error estimate for gam estimated in the same fashion as |
If a diagonal element of the Hessian is sufficiently close to 0, NA
is returned.
The cloglog.2asym
and its alias, weib.2asym
, don't converge on
occasion. This can be observed by using the trace
argument.
One strategy is to modify the initial estimates.
Kenneth Knoblauch
Klein S. A. (2001) Measuring, estimating, and understanding the psychometric function: a commentary. Percept Psychophys., 63(8), 1421–1455.
Wichmann, F. A. and Hill, N. J. (2001) The psychometric function: I.Fitting, sampling, and goodness of fit. Percept Psychophys., 63(8), 1293–1313.
glm
, optim
, glm.lambda
, mafc
#A toy example, set.seed(12161952) b <- 3 g <- 0.05 # simulated false alarm rate d <- 0.03 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- psyfun.2asym(resp.mat ~ cnt, link = probit.2asym) ddlog.glm <- psyfun.2asym(resp.mat ~ cnt, link = logit.2asym) # Can fit a Weibull function, but use log contrast as variable ddweib.glm <- psyfun.2asym(resp.mat ~ log(cnt), link = weib.2asym) ddcau.glm <- psyfun.2asym(resp.mat ~ cnt, link = cauchit.2asym) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 5) lines(pcnt, predict(ddlog.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2, col = "blue") lines(pcnt, predict(ddweib.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey") lines(pcnt, predict(ddcau.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey", lty = 2) summary(ddprob.glm)
#A toy example, set.seed(12161952) b <- 3 g <- 0.05 # simulated false alarm rate d <- 0.03 a <- 0.04 p <- c(a, b, g, d) num.tr <- 160 cnt <- 10^seq(-2, -1, length = 6) # contrast levels #simulated Weibull-Quick observer responses truep <- g + (1 - g - d) * pweibull(cnt, b, a) ny <- rbinom(length(cnt), num.tr, truep) nn <- num.tr - ny phat <- ny/(ny + nn) resp.mat <- matrix(c(ny, nn), ncol = 2) ddprob.glm <- psyfun.2asym(resp.mat ~ cnt, link = probit.2asym) ddlog.glm <- psyfun.2asym(resp.mat ~ cnt, link = logit.2asym) # Can fit a Weibull function, but use log contrast as variable ddweib.glm <- psyfun.2asym(resp.mat ~ log(cnt), link = weib.2asym) ddcau.glm <- psyfun.2asym(resp.mat ~ cnt, link = cauchit.2asym) plot(cnt, phat, log = "x", cex = 1.5, ylim = c(0, 1)) pcnt <- seq(0.01, 0.1, len = 100) lines(pcnt, predict(ddprob.glm, data.frame(cnt = pcnt), type = "response"), lwd = 5) lines(pcnt, predict(ddlog.glm, data.frame(cnt = pcnt), type = "response"), lwd = 2, lty = 2, col = "blue") lines(pcnt, predict(ddweib.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey") lines(pcnt, predict(ddcau.glm, data.frame(cnt = pcnt), type = "response"), lwd = 3, col = "grey", lty = 2) summary(ddprob.glm)
The data were obtained from the measurements of the luminance
of the R
, G
and B
channels individually,
as well as the three together, W, for each of 21 grey levels,
GL
from a screen on which a video projector was displaying
an image of a uniform field. Grey level has been normalized to
the interval [0, 1], though originally it is specified as integers
in [0, 255]. The measurements were obtained with a Photo Research
650 spectro-radiometer.
data(RGB)
data(RGB)
A data frame with 84 observations on the following 3 variables.
Lum
numeric vector of the measured luminance
in candelas/meter
GL
The grey level normalized to the interval [0, 1]
Gun
factor with levels R
G
B
W
data(RGB)
data(RGB)
This selfStart
model evaluates the parameters for describing
the luminance vs grey level relation of the R, G and B guns of
a CRT-like display, fitting a single exponent, gamma, for each
of the 3 guns. It has an initial attribute that will evaluate
initial estimates of the parameters, Blev
, Br
,
Bg
, Bb
and gamm
.
In the case of fitting data from a single gun or for a combination of guns, as in the sum of the three for calibrating the white, the parameter k
is used for the coefficient.
Both functions include gradient and hessian attributes.
SS.calib(Blev, k, gamm, GL) SS.RGBcalib(Blev, Br, Bg, Bb, gamm, Rgun, Ggun, Bgun)
SS.calib(Blev, k, gamm, GL) SS.RGBcalib(Blev, Br, Bg, Bb, gamm, Rgun, Ggun, Bgun)
Blev |
numeric. The black level is the luminance at the 0 grey level |
k |
numeric, coefficient of one gun for fitting single gun |
Br |
numeric, coefficient of the R gun |
Bg |
numeric, coefficient of the G gun |
Bb |
numeric, coefficient of the B gun |
gamm |
numeric, the exponent, gamma, applied to the grey level |
GL |
numeric, is the grey level for the gun tested, covariate in model matrix in one gun case |
Rgun |
numeric, is a covariate in the model matrix that indicates the grey level for the R gun. See the example below. |
Ggun |
numeric, is a covariate in the model matrix that indicates the grey level for the G gun |
Bgun |
numeric, is a covariate in the model matrix that indicates the grey level for the B gun |
The model
where i is in {R, G, B},
usually provides a reasonable description of the change
in luminance of a display gun with grey level, GL
.
This selfStart
function estimates and the other parameters using the
nls
function. It is assumed that grey level is normalized
to the interval [0, 1]. This results in lower correlation between
the linear coefficients of the guns, , than if the
actual bit-level is used, e.g., [0, 255], for an 8-bit graphics
card (see the example).
Also, with this normalization of the data, the coefficients,
, provide estimates of the maximum luminance for
each gun. The need for the arguments
Rgun
, Ggun
and
Bgun
is really a kludge in order to add gradient and
hessian information to the model.
returns a numeric vector giving the estimated luminance given the parameters passed as arguments and a gradient matrix and a hessian array as attributes.
Kenneth Knoblauch
~put references to the literature/web site here ~
data(RGB) #Fitting a single gun W.nls <- nls(Lum ~ SS.calib(Blev, k, gamm, GL), data = RGB, subset = (Gun == "W")) summary(W.nls) #curvature (parameter effect) is greater when GL is 0:255 Wc.nls <- nls(Lum ~ SS.calib(Blev, k, gamm, GL*255), data = RGB, subset = (Gun == "W")) MASS::rms.curv(W.nls) MASS::rms.curv(Wc.nls) pairs(profile(Wc.nls), absVal = FALSE) pairs(profile(W.nls), absVal = FALSE) #Fitting 3 guns with independent gamma's RGB0.nls <- nlme::nlsList(Lum ~ SS.calib(Blev, k, gamm, GL) | Gun, data = subset(RGB, Gun != "W")) summary(RGB0.nls) plot(nlme::intervals(RGB0.nls)) # Add covariates to data.frame for R, G and B grey levels gg <- model.matrix(~-1 + Gun/GL, RGB)[ , c(5:7)] RGB$Rgun <- gg[, 1] RGB$Ggun <- gg[, 2] RGB$Bgun <- gg[, 3] RGB.nls <- nls(Lum ~ SS.RGBcalib(Blev, Br, Bg, Bb, gamm, Rgun, Ggun, Bgun), data = RGB, subset = (Gun != "W") ) summary(RGB.nls) confint(RGB.nls)
data(RGB) #Fitting a single gun W.nls <- nls(Lum ~ SS.calib(Blev, k, gamm, GL), data = RGB, subset = (Gun == "W")) summary(W.nls) #curvature (parameter effect) is greater when GL is 0:255 Wc.nls <- nls(Lum ~ SS.calib(Blev, k, gamm, GL*255), data = RGB, subset = (Gun == "W")) MASS::rms.curv(W.nls) MASS::rms.curv(Wc.nls) pairs(profile(Wc.nls), absVal = FALSE) pairs(profile(W.nls), absVal = FALSE) #Fitting 3 guns with independent gamma's RGB0.nls <- nlme::nlsList(Lum ~ SS.calib(Blev, k, gamm, GL) | Gun, data = subset(RGB, Gun != "W")) summary(RGB0.nls) plot(nlme::intervals(RGB0.nls)) # Add covariates to data.frame for R, G and B grey levels gg <- model.matrix(~-1 + Gun/GL, RGB)[ , c(5:7)] RGB$Rgun <- gg[, 1] RGB$Ggun <- gg[, 2] RGB$Bgun <- gg[, 3] RGB.nls <- nls(Lum ~ SS.RGBcalib(Blev, Br, Bg, Bb, gamm, Rgun, Ggun, Bgun), data = RGB, subset = (Gun != "W") ) summary(RGB.nls) confint(RGB.nls)
Identical to summary.glm
but with one line of additional output: the estimate of lambda from glm.lambda
, obtained by profiling the deviance and estimating its minimum.
## S3 method for class 'lambda' summary(object, ...) ## S3 method for class 'summary.lambda' print(x, digits = max(3, getOption("digits") - 3), symbolic.cor = x$symbolic.cor, signif.stars = getOption("show.signif.stars"), ...)
## S3 method for class 'lambda' summary(object, ...) ## S3 method for class 'summary.lambda' print(x, digits = max(3, getOption("digits") - 3), symbolic.cor = x$symbolic.cor, signif.stars = getOption("show.signif.stars"), ...)
object |
Fitted model object of class “lambda” inheriting from |
x |
an object of class “summary.lambda”, usually a result of a call to |
digits |
the number of significant digits to use when printing. |
symbolic.cor |
logical. If |
signif.stars |
logical. If |
... |
further arguments passed to or from other methods. |
Provides a summary of the class lambda
object generated by glm.lambda
.
Returns the same structure as summary.glm
with an added component, lambda
. is the estimated upper asymptote of the psychometric function.
Ken Knoblauch