Understanding the Jeffrey prior penalty term in the glmfit function in matlab R2024a

Hi,

Consider the example separable dataset

x_i = 800:5:900;
y_i = [0 0 0 0 0, ...
    0 0 0 0 0, ...
    1 1 1 1 1, ...
    1 1 1 1 1, ...
    1];

In the R2024a version of matlab it is now possible to fit a probit curve to these data using Jeffrey's prior as a penalty function.

%Standard way - singular coefficients
b_standard = glmfit(x_i', y_i', 'binomial','link', 'probit');

Warning: Iteration limit reached.

Warning: The estimated coefficients perfectly separate failures from successes. This means the theoretical best estimates are not finite. For the fitted linear combination XB of the predictors, the sample proportions P of Y=N in the data satisfy:
XB<-0.0197012: P=0
XB>-0.0197012: P=1
To get proper estimates, fit with 'LikelihoodPenalty' set to 'jeffreys-prior'.

%Jeffrey prior - nonsingular coefficients
b_jeff = glmfit(x_i', y_i', 'binomial','link', 'probit',LikelihoodPenalty="jeffreys-prior");

%Plot the standard and Jeffrey-prior coefficients
figure(1)
hold on
plot(b_standard(1),b_standard(2),'r*')
plot(b_jeff(1),b_jeff(2),'k*')

This is good because instead of obtaining divergent coefficients we obtain finite coefficients which are shifted towards the origin.

My questions are related to what is going on under the hood of the glmfit function.

• Is there a simple way to plot the Jeffrey prior that is being used (as a function of possible coefficient values), and if so how?
• Is it also possible to plot the resulting likelihood function (as a function of possible coefficient values) when we used the Jeffrey prior?

In the case where the Jeffrey prior is not being used, I can compute the likelihood function as

%Plot the likelihood
%I here transformed from b_0 and b_1 to the mean and slope parameters
%because the plots often look nicer in the latter coordinates
mu = linspace(840,860,201);
sigma = linspace(-1,5,101);

LikelihoodList = [];
for i = 1:length(mu)
        for j = 1:length(sigma)
            LikelihoodList = [LikelihoodList; mu(i), sigma(j), LikelihoodFunction(-mu(i)./sigma(j),1./sigma(j),x_i,y_i)];
        end
end
X = LikelihoodList(:,1); Y = LikelihoodList(:,2); [x,y]=meshgrid(X,Y); Z = LikelihoodList(:,3);
figure(2)
plot3(X,Y,Z)

function likelihood = LikelihoodFunction(b_0,b_1,x_i,y_i)
    %Computes the likelihood
    %Take in a single value of b0 and b1 and a list of experimental data x_i and y_i
    p_i = normcdf(b_0 + b_1*x_i);
    likelihood_product_list = p_i.^y_i.*(1-p_i).^(1-y_i);
    likelihood = prod(likelihood_product_list);
end

Related texts can be found here:

brglm: Bias Reduction in Binomial-Response Generalized Linear Models (r-project.org)

brglm2: Bias Reduction in Generalized Linear Models (r-project.org)

[P] Fit Logistic Regression with Jeffreys Prior : r/MachineLearning (reddit.com)

Firth's bias correction as a Bayesian model | Random effect

neither are using matlab.