Hi Craig,
To compute the confidence limits about the quantile ( q ) for a logistic regression model, you need to invert the logistic function and then calculate the confidence interval for the quantile. Here’s a step-by-step approach to achieve this in MATLAB:
- Fit the logistic regression model.
- Calculate the confidence intervals for the model coefficients.
- Determine the quantile confidence intervals by inverting the logistic function using the confidence intervals of the coefficients.
Below is the MATLAB code to achieve this:
% Fit the logistic regression model
xWT = [5.5,16.5,11,13.8,10.1,14.7,10.4,11.7,9.7,7.3,7.8,8.1,12.2,8.5,11.8,11.7121,11.4083,11.1558,12.4633,12.2761,12.1107,11.9628,11.8291,11.7072,11.5952,11.4917,11.3955,11.3057,11.2214,11.1421]';
yWT = [0,1,0,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1]';
n = ones(size(yWT));
[b, dev, stats] = glmfit(xWT, [yWT, n], 'binomial', 'Link', 'probit');
% Define the range for prediction
wt = 5 : 0.05 : 17;
% Calculate the fitted values and their confidence intervals
[y_fit, dylo, dyhi] = glmval(b, wt, 'probit', stats, 'confidence', 0.95);
% Calculate the confidence intervals for the coefficients
ci = coefCI(fitglm(xWT, yWT, 'Distribution', 'binomial', 'Link', 'probit'), 0.95);
% Invert the logistic function to find the quantiles
% For a given probability p, the quantile q is given by:
% q = -b(1) / b(2) + (log(p / (1 - p))) / b(2)
% Define the range of probabilities
p = 0.01:0.01:0.99;
% Calculate the quantiles
q = -b(1) / b(2) + norminv(p) / b(2);
% Calculate the confidence intervals for the quantiles
q_lo = -ci(1,2) / ci(2,2) + norminv(p) / ci(2,2);
q_hi = -ci(1,1) / ci(2,1) + norminv(p) / ci(2,1);
% Plot the results
figure;
hold on;
plot(xWT(yWT == 1), yWT(yWT == 1), 'ko', 'LineWidth', 0.5, 'MarkerFaceColor', 'g', 'MarkerSize', 8);
plot(xWT(yWT == 0), yWT(yWT == 0), 'ks', 'LineWidth', 0.5, 'MarkerFaceColor', 'r', 'MarkerSize', 8);
plot(wt, y_fit, '-m');
plot(wt, y_fit - dylo, '-k');
plot(wt, y_fit + dyhi, '-k');
plot(q, p, '--r');
plot(q_lo, p, '--r');
plot(q_hi, p, '--r');
xlim([5, 17]);
title('Wu and Tian Example');
xlabel('Quantile (q, in units wt)');
ylabel('Probability (p)');
grid on;
legend({'Go', 'No-Go', 'Logistic Fit', 'Confidence Bounds from glmval', 'Quantile Confidence Interval'}, 'Location', 'east');
Explanation:
- The glmfit function fits a logistic regression model to the data.
- The glmval function computes the predicted values and their confidence intervals.
- The coefCI function calculates the confidence intervals for the model coefficients.
- The quantiles are computed by inverting the logistic function using the model coefficients.
- The confidence intervals for the quantiles are calculated using the confidence intervals of the model coefficients.
In the plot, the magenta line represents the logistic fit, the black lines are the confidence bounds from glmval, and the red dashed lines are the quantile confidence intervals. This should match the results given by Gonogo.
