# Model-free estimation of a psychometric function

## Example 1. Frequency of seeing

Miranda, M. A. & Henson, D. B. “Perimetric sensitivity and response variability in glaucoma with single-stimulus automated perimetry and multiple-stimulus perimetry with verbal feedback”, Acta Ophthalmologica, 86, 202-206, 2008.

Back to Examples

A flash of light of variable intensity was presented repeatedly at a fixed location in the visual field of a subject who reported whether the flash was visible. There were 3–20 trials at each stimulus level.

 `clear, load examples/example_01;`

Loaded MatLab `data` consist of three columns:

1. Stimulus level, `x`
2. Number of successes, `r`
3. Number of trials, `m`
Parametric fitting

Four different parametric models are fitted to these data: Gaussian (probit), Weibull, reverse Weibull and logistic. The parameters of the models are estimated in the Matlab functions `binomfit_lims` for the probit and logit links. The other two models require an estimate of the exponent `K` and the estimation is performed in the functions `binom_weib` and `binom_revweib`. The values of the fitted functions at specified points are calculated in the function `binomval_lims`.

First plot the psychometric data (black dots):

 `figure; plot( x, r ./ m, 'k.'); axis([0.05 1.35 -0.02 1]); axis square;`

1. For the Gaussian cumulative distribution function (black curve):

 ```degpol = 1; % Degree of the polynomialb = binomfit_lims( r, m, x, degpol, 'probit' ); numxfit = 999; % Number of points to be generated minus 1 xfit = [min( x ):(max(x)-min(x))/numxfit:max( x )]'; % Plot the fitted curve pfit = binomval_lims( b, xfit, 'probit' ); hold on, plot( xfit, pfit, 'k' );```

2. For the Weibull function (red curve):

 ```[ b, K ] = binom_weib( r, m, x ); guessing = 0; % guessing rate lapsing = 0; % lapsing rate % Plot the fitted curve pfit = binomval_lims( b, xfit, 'weibull', guessing, lapsing, K ); hold on, plot( xfit, pfit, 'r' );```

3. For the reverse Weibull function (green curve):

 ```[ b, K ] = binom_revweib( r, m, x ); % Plot the fitted curve pfit = binomval_lims( b, xfit, 'revweibull', guessing, lapsing, K ); hold on, plot( xfit, pfit, 'g' );```

4. For the logistic function (blue curve):

 ```b = binomfit_lims( r, m, x, degpol, 'logit' ); % Plot the fitted curve pfit = binomval_lims( b, xfit, 'logit' ); hold on, plot( xfit, pfit, 'b' );```

Local linear fitting

Local linear fitting is performed in the function `locglmfit`, which returns the fitted values at specified points. This function requires a bandwidth as an input. The bandwidth `bwd` is typically chosen by a cross-validation. There are three different loss functions used in cross-validation: ISE defined on a p-scale, ISE defined on an eta-scale, and deviance.

 ```bwd_min = min( diff( x ) );bwd_max = max( x ) - min( x );bwd = bandwidth_cross_validation( r, m, x, [ bwd_min, bwd_max ] );```

The values of cross-validation bandwidths are

• `bwd = 0.1000` for ISE on a p-scale,
• `bwd = 0.3128` for ISE defined on an eta-scale,
• `bwd = 0.2959` for deviance.

Here, the bandwidth obtained with cross-validated deviance is used.

 ```bwd = bwd(3); % Choose the third estimate, which is based on cross-validated deviance pfit = locglmfit( xfit, r, m, x, bwd ); figure; plot( x, r ./ m, 'k.'); axis([0.05 1.35 -0.02 1]); axis square; hold on, plot( xfit, pfit, 'k' );```

From the fitted values of the psychometric function, the threshold and slope for a requied threshold level, here `prob` = 0.5, are calculated in the function `threshold_slope`. The functions `bootstrap_sd_th` and `bootstrap_sd_sl` estimate their standard deviations by the percentile bootstrap method. In the following example, 200 bootstrap replications were used. Note that these functions, as with all bootstrap methods, return similar but slightly different values at each execution.

 ```prob = 0.5; % Required threshold level niter = 200; % Number of bootstrap iterations [ threshold, slope ] = threshold_slope( pfit, xfit, prob ); sd_th = bootstrap_sd_th( prob, r, m, x, niter, bwd ); % Be patient, slow process sd_sl = bootstrap_sd_sl( prob, r, m, x, niter, bwd ); % Be patient, slow process```

Examples of values returned by MatLab are as follows:

• `threshold (sd_th) = 0.9745 (0.0558)`
• `slope (sd_sl) = 1.5536 (0.2941)`

Bootstrap estimates of confidence intervals for the threshold and slope are calculated by the functions `bootstrap_ci_th` and `bootstrap_ci_sl;` here a significance level `alpha` = 0.05 was used. Again, the these functions return similar but slightly different values at each execution.

 ```prob = 0.5; % Required threshold level alpha = 0.05; % Significance level for the confidence intervals niterci = 1000; % Number of bootstrap iterations ci_th = bootstrap_ci_th( prob, r, m, x, niterci, bwd, alpha); % Be patient, slow process ci_sl = bootstrap_ci_sl( prob, r, m, x, niterci, bwd, alpha ); % Be patient, slow process```

Examples of confidence intervals returned by MatLab are as follows:

• `ci_th = [0.8363,1.0538];`
• `ci_sl = [0.9863,2.1536];`
Back to Examples