How to test if the difference between two logit-function-fits is significant or not?

Hi,

I am considering a binomial experiment where each trial can only have two possible outputs (0,1). I vary the value of a predictor (experiment parameter) and repeat the experiment several times for each predictor value. Averaging the responses over trials for each predictor value to obtain the probability of outcome=1 as a function of the predictors, yields s-shaped curves like the blue and red curves below (depending on the experiment fixed conditions):

I would like to test if the two curves are significantly different and quantify the bias between them.

I tried to solve this issue in my way and it works, I would be very grateful if anybody could confirm this is the good way to do it or there is a better one. Furthermore I would really appreciate any advice on how to calculate the optimal sample size for this kind of analysis. First I will try to explain in a general way what I did and then I will show my MATLAB implementation.

General idea

I assume the relationship between a vector of predictors x = [x1,x2,..xi,xn] and the vector of success probabilities P=[P(Y=1|xi), 1<=i<=n ] could be modelled by logit link function as:

P=1/ 1+exp[(b1+ (b2*x)]; (1)

where b1 and b2 determine the logistic intercept and slope.

To test if the difference in intercept or slope btw the two datasets is significant I will add two other parameters to the equation {b3,b4} that will depend on a constant 'I' that will be I=0 for the fit of the first dataset and I=1 to fit the second one

P=1/ 1+exp[(b1+ (b2*x) + (b3*I) + (b4*x*I)]; (2)

Therefore I will fit the first dataset obtaining a maximum likelihood estimation (MLE) of b1 and b2 and their corresponding p-values using the model in equation (1) and assuming inter-trials binomial variability. Then I will fit the second dataset with the model (2) by fixing b1 and b2 as constants (obtained from the first estimation) to get an estimate of b3 and b4 and their p-values.

If b3 and b4 are significantally different from zero (pvalue<0.05), then the two datasets yielding to the blue and red curves are significantilly different.

If the p-value of b3 is <0.05 and the p-value of b4 is >0.05 The difference is the result of an horizontal offset. Therefore I can use eq.2 to fit again the second dataset this time removing b4 (eq3):

P=1/ 1+exp[(b1_estimate+ (b2_estimate*x) + (b3*I)]; (3)

Now b3 relate to the horizontal offset. To quantify this offset in predictor (x) units I have to divide it by the slope parameter

offset = b3/b2 ; (4)

Matlab implementation

%First I generate two sets of synthetic data for the test (y1 & y2)

x = -50:5:50; % predictors
n = 20; % num of samples for each predictor

% Dataset 1
beta1 = [3 ,0.2]; % parameters logit function,...
% beta(1)/beta(2) = intercept; beta(2) = ""slope"";  
[yfit1 ] = glmval(beta1',x ,'logit'); % true generative distribution
y1  = binornd(n,yfit1); % binomial sampling from the true distribution

% Dataset 2
beta2 = [0.5 ,0.2]; 
[yfit2 ] = glmval(beta2' ,x,'logit'); 
y2  = binornd(n,yfit2); 

%Display simulated data for the two curves
cla
plot(x, y1./n, 'ob', x, yfit1, '--b','linewidth',2)
hold on
plot(x, y2./n, 'or', x, yfit2, '--r','linewidth',2)
set(gca,'xtick',[-50:10:50])
ylabel('P')
xlabel('Predictors')

%  Dataset1 fit  
 [b, dev1, stat1] = glmfit(x, [y1 repmat(n,size(y1)) ],...
'binomial','link', 'logit');

   % Dataset2 fit  - to add the offset corresponds to fix b1 and b2 and determine
   % b3 and b4
   W = b(1)+(b(2)*x); 
   [tmp, dev2, stat2] = glmfit(x, [y2 repmat(n,size(y2))],...
 'binomial','link','logit','offset',W','constant','on');
   b([3,4]) = tmp; 

   %p-value of b(3)  and b(4)
   b3_pvalue = stat2.p(1)
   b4_pvalue = stat2.p(2)

   % if for b(3) p<0.05 and for b(4) p>0.05 I remove b(4) from the equation
    [offsetFit dev3  stat3]=glmfit(ones(size(y2)),...
[y2 repmat(n,size(y2))],'binomial','link','logit',...
'const','off','offset',W') ;

    % offset in predictor unit
    offset = offsetFit/b(2)