glm with probit model

[Daniel]

I have successfully implemented a glmfit with binary response and a probit model. My model has the form of p(t)=probit(-x^2). This seems to work fairly well until a lapse occurs and subjects miss a level that they "should" have detected. For example, maybe they fall asleep for a few moments. I understand that this can be modeled by including a lapse rate lambda (usually assume to be less than 0.05 or so) yielding a model of the form:

lambda + (1 - lambda)*probit (-x^2)

I haven't been able to figure out how to fit this with glmfit. Am I missing something? Ideas?

[Peter Perkins]

Daniel, your description is a little vague to me. It sounds like you're using the GLMFIT function in the Statistics Toolbox with a probit link and a binomial distribution. It sounds like what you want is to model the probability of "success" (whatever that means in your case) as

p = lambda + (1-lambda)*Phi(eta)

where Phi is the standard normal CDF, and eta is the linear predictor X*beta. Is that right?

You can do this by defining your own link function as described in the help for GLMFIT. I believe the usual thing described in the literature is to fix lambda, fit the GLM and compute the deviance, and then vary lambda to find the value that minimizes the deviance.

The one catch in doing this with GLMFIT is that you'll probably find that the standard GLM initialization doesn't work well with that kind of link function -- depending on your data and the value of lambda, it may try to initialize with probabilities that are less than lambda. That's easy to fix if you can modify GLMFIT in your MATLAB installation.

If you can't, there was a thread a few weeks ago with the subject "glmfit binomial-probit with natural mortality" that describes a way to call GLMFIT to fit this model.

Hope this helps.

- Peter Perkins
The MathWorks, Inc.

[Daniel]

Peter, thanks for your reply. Yes, you seem to have understood my model even though my description was not clear. I understand that I can create my own functions, and I understand the example that you provided. However, I am still struggling with how to implement my model, so let me first spell out my model more clearly. The underlying model is based on a probit link. To make this even more explicit, here is the call that I use for model fitting:
[b,dev2,stats]=glmfit(motion,correct,'binomial','link','probit');

In fact, this works very well for much of my data. However, humans sometimes "miss" on trials when the probability says that they should almost never miss and this can severely bias the desired parameter fits from the model. (An example of a lapse is that I should hear that tone except that I fell asleep while you played the tone for me.) This is often handled in the psychophysical literature by adding one additional parameter to the model to represent the fact that subjects might sometimes miss on trials that they should have detected. (See below for an equation showing how this is implemented.)

To avoid ambiguity, I'll provide some portions of an m-file that will explicitly outline my model.

mu=-3; %This is the location of the peak of the probability distribution function.
sigma=3; %This is the standard deviation.
lambda=0.01; % This is my lapse rate.

x=-20:.01:20;

If I understand the m-file you provided, the following represent the three functions that are required.
deta=(1-2*lambda)*normpdf(x,mu,sigma); %This is the derivative of the link function
eta=norminv((p_lapse-lambda)/(1-2*lambda),mu,sigma); %This is a form of the link function
mu=lambda + (1-2*lambda)*normcdf(x,mu,sigma); %This is a form of the inverse link function

However, when I try to use the m-file that you provided as a template, I start to run into problems.

I can almost figure out how to implement the derivative function, since I can write this equation in a simple closed form solution:
deta=(1-2*L)*exp(-((mu)).^2)/sqrt(2*pi);
(But this is missing the normalization by sigma which I can't see how to write in the desired format for the GLM, but I probably am just missing something.)

However, the inverse link function cannot be written in a closed form solution. In fact it includes the normal cdf, which is a from of the error function and I just can't see how to write this in a closed form that is appropriate for the glmfit. Seems like the same is true for the link function.

Am I missing something?

Thanks in advance for any advice. Sincerely,
Dan

[Peter Perkins]

Daniel wrote:

> If I understand the m-file you provided, the following represent the three functions that are required.
> deta=(1-2*lambda)*normpdf(x,mu,sigma); %This is the derivative of the link function
> eta=norminv((p_lapse-lambda)/(1-2*lambda),mu,sigma); %This is a form of the link function
> mu=lambda + (1-2*lambda)*normcdf(x,mu,sigma); %This is a form of the inverse link function
>
> However, when I try to use the m-file that you provided as a template, I start to run into problems.
>
> I can almost figure out how to implement the derivative function, since I can write this equation in a simple closed form solution:
> deta=(1-2*L)*exp(-((mu)).^2)/sqrt(2*pi);
> (But this is missing the normalization by sigma which I can't see how to write in the desired format for the GLM, but I probably am just missing something.)
>
> However, the inverse link function cannot be written in a closed form solution. In fact it includes the normal cdf, which is a from of the error function and I just can't see how to write this in a closed form that is appropriate for the glmfit. Seems like the same is true for the link function.

I'm not sure where those 2's are coming from.

These are the probit link, link derivative, and inverse link functions that GLMFIT uses (edit private/stattestlink, called by glmfit):

link = @(mu) norminv(mu);
dlink = @(mu) 1 ./ normpdf(norminv(mu));
% keep mu = ilink(eta) in [eps, (1-eps)];
lowerBnd = norminv(eps(dataClass)); upperBnd = -lowerBnd;
ilink = @(eta) normcdf(min(max(eta,lowerBnd),upperBnd));

Given lambda, I believe you want

ilink = @(eta) lambda + ...
(1-lambda).*normcdf(min(max(eta,lowerBnd),upperBnd));

and therefore

link = @(mu) norminv((mu-lambda)./(1-lambda));
dlink = @(mu) 1 ./ ((1-lambda).*normpdf(norminv((mu-lambda)./(1-lambda))));

but you'll want to check my work. If your lambda is really .01, then you may not run into the problem with initial values that I alluded to, unless youhave more than 100 observations.

Hope this helps.

[Daniel]

Thanks very much. We will try your suggestions and let you know.

FYI, the 2 comes from the fact that people can lapse on both ends, so we want the cdf to range from lambda on the low end to (1-lambda) on the high end. For example, if lambda equals 0.01, then the lowest probability for very small values of the independent variable would be 0.01 and the highest probability for very large values of the independent variable would be 0.99.

Hope that makes sense. Will let you know if this works.

[Peter Perkins]

Daniel wrote:

> FYI, the 2 comes from the fact that people can lapse on both ends, so we want the cdf to range from lambda on the low end to (1-lambda) on the high end. For example, if lambda equals 0.01, then the lowest probability for very small values of the independent variable would be 0.01 and the highest probability for very large values of the independent variable would be 0.99.
>
> Hope that makes sense. Will let you know if this works.

It does, but my versions of those three functions are not what you want. But I think the changes needed are obvious.

[Daniel]

Hi Peter , thanks for all of your help. it seems that I am very close but I seem to missing something critical. Following your suggestion, I defined the following functions:

myilink = @(eta) L + (1-2*L).*normcdf(min(max(eta,lowerBnd),upperBnd));
mylink = @(mu) norminv((mu-L)./(1-2*L));
mydlink = @(mu) 1 ./ ((1-2*L).*normpdf(norminv((mu-L)./(1-L))));

If I fit with the standard probit model:
b = glmfit(motion,y,'binomial','link','probit');
everything works fine, so I know that my input vectors and general approach are correct.

However, when I fit using the following call (as I believe you had suggested):
b = glmfit(motion,y,'binomial','link',{mylink mydlink myilink});
I get the following error message:

??? Input argument "L" is undefined.

Error in ==> direction_detect>@(mu,L)norminv((mu-L)./(1-2*L)) at 64
mylink = @(mu,L) norminv((mu-L)./(1-2*L));

Error in ==> glmfit at 277
eta = linkFun(mu);

Error in ==> direction_detect at 77
b = glmfit(motion,y,'binomial','link',{mylink mydlink myilink});

First, I should probably let you know that L is defined and is a scalar. It is defined via the following three lines of m-file code:
lambda=0:.01:.06;
for j=1:length(lambda)
L = lambda(j);
so I believe that it both exists and is a scalar.

Now, if were to be calling my function(s) directly myself, I know that I could simply add my variable L in the function definition. For example:
myilink = @(eta, L) L + (1-2*L).*normcdf(min(max(eta,lowerBnd),upperBnd));

However, since I do not control the call function within glmfit, I do not see how to pass the L variable to my function(s).

Suggestions?

Thanks in advance,
Dan

[Tom Lane]

> myilink = @(eta) L + (1-2*L).*normcdf(min(max(eta,lowerBnd),upperBnd));
> mylink = @(mu) norminv((mu-L)./(1-2*L));
> mydlink = @(mu) 1 ./ ((1-2*L).*normpdf(norminv((mu-L)./(1-L))));

...

> I get the following error message:
>
> ??? Input argument "L" is undefined.

Dan, you should set L before you define the functions above. Its value is
then incorporated into the function definitions.

It sounds like you're setting L in a loop. Then re-define these functions
each time through the loop, so they pick up the new value of L. I don't
think there will be a lot of overhead in doing that.

-- Tom

[Peter Perkins]

Daniel wrote:
> Hi Peter , thanks for all of your help. it seems that I am very close but I seem to missing something critical. Following your suggestion, I defined the following functions:
>
> myilink = @(eta) L + (1-2*L).*normcdf(min(max(eta,lowerBnd),upperBnd));
> mylink = @(mu) norminv((mu-L)./(1-2*L));
> mydlink = @(mu) 1 ./ ((1-2*L).*normpdf(norminv((mu-L)./(1-L))));

I think you want 1-2*L in the denom of the last one.

> However, when I fit using the following call (as I believe you had suggested):
> b = glmfit(motion,y,'binomial','link',{mylink mydlink myilink});
> I get the following error message:
>
> ??? Input argument "L" is undefined.
>
> Error in ==> direction_detect>@(mu,L)norminv((mu-L)./(1-2*L)) at 64
> mylink = @(mu,L) norminv((mu-L)./(1-2*L));

Where did that link function come from? It isn't the one you cited above, and it has L as an argument. L should not be an input argument. Let the anonomous function "captures" L by creating a new mylink at each iteration inside this loop:

[Daniel]

Peter Perkins <Peter.Perki...@mathworks.com> wrote in message <gfcbo7$sfi$1...@fred.mathworks.com>...

> > mydlink = @(mu) 1 ./ ((1-2*L).*normpdf(norminv((mu-L)./(1-L))));
>
> I think you want 1-2*L in the denom of the last one.
>

You are correct. My mistake.

> > mylink = @(mu,L) norminv((mu-L)./(1-2*L));
>
> Where did that link function come from? It isn't the one you cited above, and it has L as an argument. L should not be an input argument. Let the anonomous function "captures" L by creating a new mylink at each iteration inside this loop:
>
> > lambda=0:.01:.06;
> > for j=1:length(lambda)
> > L = lambda(j);

You are correct. I tried using the function above after I had tried:

mylink = @(mu) norminv((mu-L)./(1-2*L));

Both gave the same error message.

I think that it is working now. When I am certain, would you like the "final" code. I simply modified the example code that you had provided earlier with the updated model and model fit, etc.

[Peter Perkins]

Daniel wrote:

> I think that it is working now. When I am certain, would you like the "final" code. I simply modified the example code that you had provided earlier with the updated model and model fit, etc.

Post it on the MATLAB Central File Exchange!

[Matthew Crema]

Hi Daniel and Peter,

I realize this is an old thread, but was the "final" code for this problem posted somewhere?

As in the problem worked above, it is also common in experimental psychology to fit a "psychometric function" of the same form to binomial data:
G + (1-G-L).*normcdf(x);

where G is the "Guess Rate" and L is the "Lapse Rate". (In your example G=L)

After reading this thread, I think I understand how to perform the appropriate fit in MATLAB, but I get beta = [NaN; NaN] when I run my code (below) using non-zero G and L values.

If I set G and L to zero, the code returns a finite value for beta (the same beta I would get with link='probit') but the goodness of fit is of course poor.

Thank you.

-----
% Mock data generated using G=1/3, L=0, mu = 60, sigma = 10;
Gact = 1/3;
Lact = .001;
mu = 60;
sigma = 10;
X = [24;31;31;38;39;43;43;48;49;56;67;67;70;73;73;77;77;78;87;88];
Y = [3;2;2;3;1;2;4;0;4;1;6;5;6;4;7;9;1;1;10;7];
N = [9;6;5;9;4;5;10;1;10;2;7;6;7;4;7;9;1;1;10;7];

% Use G = 0 and L = 0 for glmfit (this is wrong, want G = Gact, L = Lact)
G = 0;
L = 0;

% Link function (inv normal transformed to approach guess and lapse rates)
link = @(mu) norminv((mu-G)./(1-G-L));
% Derivative of link function
dlink = @(mu) 1 ./ ((1-G-L).*normpdf(norminv((mu-G)./(1-G-L))));
lowerBnd = norminv(eps('double')); upperBnd = -lowerBnd;
% Inverse link
ilink = @(eta) G + (1-G-L).*normcdf(min(max(eta, lowerBnd), upperBnd));

% Perform GLM fit
[beta, ~, stats] = glmfit(X, [Y N], 'binomial', 'link', {link, dlink, ilink});

% Plot Results
plot(X, Y./N, 'o', X, glmval(beta, X, {link, dlink, ilink}, 'size', N)./N, ...
X, Gact + (1-Gact-Lact)*normcdf(X, mu, sigma), 'k--')
legend('Raw Data', 'GLM Fit', 'Actual PMF', 'Location', 'SouthEast')
-----

[Matthew Crema]

"Matthew Crema" wrote in message <lnpqns$4ed$1...@newscl01ah.mathworks.com>...

> Hi Daniel and Peter,
>
> I realize this is an old thread, but was the "final" code for this problem posted somewhere?
>

I think I have a solution after reviewing another thread on the subject "glmfit binomial-probit with natural mortality" There is some code that assumes the curve approaches an unknown baseline value and determines the best-fit curve by maximum likelihood estimation.

My "final" code, below, assumes the asymptotes are known and forces the curve to approach these.

Best,
-Matt

-----------
%% MYGLMFIT Fit GLM to Binomial Data with non-zero Guess and Lapse Rates
% B = MYGLMFIT(X, Y, N, FITTYPE, G, L) fits a generalized linear
% model to binomial data using the predictor matrix X, given the number
% of hits Y and total number of observations N for each entry in X. B is
% a vector of coefficient estimates.
%
% FITTYPE must be 'probit' or 'logit'. Unlike the builtin GLMFIT
% function, MYGLMFIT assumes that the Guess Rate G (for small X), and
% Lapse Rate L (for large X) are known and are generally non-zero.
%
% [B, YHAT, DYLO, DYHI] = YGLMFIT(X, Y, N, FITTYPE, G, L) computes
% predicted values for the model YHAT and computes 95% confidence bounds
% on the predicted Y values. DYLO and DYHI define a lower confidence
% bound of YHAT-DYLO and an upper confidence bound of YHAT+DYHI.
%
% See also GLMFIT, GLMVAL.
%
% Written by Matthew Crema - 2014 borrowing some code and ideas from Tom
% Lane and Peter Perkins found on MATLAB newsgroup.

function [beta_hat, yhat, dylo, dyhi] = myglmfit(X, Y, N, fittype, G, L)

%% Check input arguments

% Determine Guess and Lapse Rates
if ~exist('G', 'var')
G = 0; % Assume Guess Rate is 0
end
if ~exist('L', 'var')
L = 0; % Assume Lapse Rate is 0
end

% Deal with situation where Y passed as two column matrix
if isequal(size(Y,2), 2) && ~exist('N', 'var')
N = Y(:,2); % Number of observations
Y = Y(:,1); % Number of hits
end

% Determine which link function to use
if ~exist('fittype', 'var') || isequal(fittype, 'probit')
mylink = {@probitlink, @probitdlink, @probitilink};
elseif isequal(fittype, 'logit')
mylink = {@logitlink, @logitdlink, @logitilink};
else
error('MYGLMFIT:FITTYPE', 'Unrecognized Fit Type')
end

%% Compute GLM Fit

% Set INIT flag to indicate first iteration, so that when GLMFIT is called,
% Link functions will compute an initial guess.
init = true;

if isequal(nargout, 1)
% Parameter estimation only
beta_hat = glmfit(X, [Y N], 'binomial', 'link', mylink);
else
% Estimate parameters and return best fit and confidence limits
[beta_hat, ~, stats] = glmfit(X, [Y N], 'binomial', 'link', mylink);

if ~all([isnan(stats.se(:)); isnan(stats.coeffcorr(:))])
% Get YHAT and confidence intervals
[yhat, dylo, dyhi] = glmval(beta_hat, X, mylink, stats);
else
% Can not compute confidence intervals. Avoid an error in GLMVAL.
yhat = glmval(beta_hat, X, mylink, stats);
dylo = nan(size(yhat));
dyhi = nan(size(yhat));
end
end

%% Nested Functions Defining Links for Logit and Probit Fits

% Logit Link Function with non-zero Guess and Lapse Rates
function eta = logitlink(mu)
% Link function
if init
% Adjust starting values first iteration
mu = min(max(mu, G + (1 - G)*.5/2), 1 - (1 - G)*.5/2);
end
% Standard Logit link
%eta = log(mu ./ (1-mu));

% Logit link with non-zero Guess and Lapse Rates
eta = log((mu-G)./(1-L-mu));
end

function deta = logitdlink(mu)

% Derivative of link function

if init
% Adjust starting values first iteration
mu = min(max(mu, G + (1 - G)*.5/2), 1 - (1 - G)*.5/2);
init = false;
end
% Standard dLogit link
%deta = 1 ./ (mu .* (1-mu));

% dLogit link with non-zero Guess and Lapse Rates
deta = (1-G-L)./((mu-G).*(1-L-mu));
end

function mu = logitilink(eta, lowerBnd, upperBnd)
% Inverse of link function
if nargin < 2
% keep mu = ilink(eta) in (approx) [G+eps, (1-eps)];
lowerBnd = log(eps);
upperBnd = -lowerBnd;
end

% Standard iLogit link
%mu = 1./(1 + exp(-min(max(eta, lowerBnd), upperBnd)));

% iLogit link with non-zero Guess and Lapse Rates
mu = G + (1-G-L)./(1 + exp(-min(max(eta, lowerBnd),upperBnd)));
end

% Probit Link Function with non-zero Guess and Lapse Rates
function eta = probitlink(mu)
% Link function
if init
% Adjust starting values first iteration
mu = min(max(mu, G + L + .01), 1 - L);
end
% Standard Probit link
%eta = norminv(mu);

% Probit link with non-zero Guess and Lapse Rates
eta = norminv((mu-G)./(1-G-L));
end

function deta = probitdlink(mu)

% Derivative of link function

if init
% Adjust starting values first iteration
mu = min(max(mu, G + L + .01), 1 - L);
init = false;
end
% Standard dProbit link
%deta = 1./normpdf(norminv(mu));

% dProbit link with non-zero Guess and Lapse Rates
deta = 1./((1-G-L).*normpdf(norminv((mu-G)./(1-G-L))));
end

function mu = probitilink(eta, lowerBnd, upperBnd)
% Inverse of link function
if nargin < 2
% keep mu = ilink(eta) in (approx) [G+eps, (1-eps)];
lowerBnd = norminv(eps);
upperBnd = -lowerBnd;
end

% Standard iProbit link
%mu = 1./(1 + exp(-min(max(eta, lowerBnd), upperBnd)));

% iProbit link with non-zero Guess and Lapse Rates
mu = G + (1-G-L).*normcdf(min(max(eta, lowerBnd), upperBnd));
end
end


-----------