Path: news.mathworks.com!not-for-mail
From: "Daniel " <dan.merfeld@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: glm with probit model
Date: Sun, 9 Nov 2008 09:58:01 +0000 (UTC)
Organization: Harvard University
Lines: 31
Message-ID: <gf6c79$87j$1@fred.mathworks.com>
References: <geapba$62l$1@fred.mathworks.com> <gecis6$hje$1@fred.mathworks.com>
Reply-To: "Daniel " <dan.merfeld@gmail.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1226224681 8435 172.30.248.35 (9 Nov 2008 09:58:01 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 9 Nov 2008 09:58:01 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 45502
Xref: news.mathworks.com comp.soft-sys.matlab:499731


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