Effacer les filtres
Effacer les filtres

Compute Confidence Interval about q of logistic binary fit

3 vues (au cours des 30 derniers jours)
Craig
Craig le 26 Avr 2024
Commenté : Craig le 26 Juin 2024
I have results from a binary study. The study was analyzied by others using Gonogo (https://www2.isye.gatech.edu/~jeffwu/sensitivitytesting/), the prefered code developed by/for the government for such studies. My code for anlaysis is shown below. I am able to fit a logistic curve using glmfit (magentia line in plot). I then compute the confidence limits about the probablility p using glmval (black lines). These black line match well with the results given by Gonogo (dashed blue lines). However, I am unable to figure out how to compute the confidence limits about the quantile q (the Gonongo results are shown as dashed red lines). I tried using coefCI, with the results in cyan. It is doing something but I don't understand exactly what, and it does not really match the dashed red lines that I am looking for. Does anyone know how I can compute the confidence limits about the quantile q? Your assistance is appreciated.
%% Output from Gonogo - Main code in next section
R_mu = 10.1707909;
R_sigma = 0.9344091;
R_b0 = -10.88473;
R_b1 = 1.070195;
R_output = [
1.124021 5.5 9.875979 0 0 0.376189
1.717474 5.833333 9.949193 0 2.0e-06 0.406269
2.310404 6.166667 10.022929 0 9.0e-06 0.437133
2.902698 6.5 10.097302 0 4.3e-05 0.468657
3.494207 6.833333 10.17246 0 0.000177 0.500713
4.084734 7.166667 10.248599 0 0.000652 0.533182
4.477745 7.388889 10.300033 0 0.001455 0.555004
5.065991 7.722222 10.378454 0 0.004391 0.587936
5.457105 7.944444 10.431784 0 0.008595 0.609998
5.847164 8.166667 10.486169 2.0e-06 0.015984 0.632136
6.235924 8.388889 10.541854 1.3e-05 0.028261 0.654357
6.490555 8.534934 10.57931 4.1e-05 0.040 0.669016
6.808555 8.717996 10.62744 0.00016 0.060 0.687473
7.008118 8.833333 10.658548 0.000356 0.076166 0.699163
7.199693 8.944444 10.689196 0.000737 0.094688 0.710483
7.390505 9.055556 10.720606 0.001463 0.116333 0.721871
7.571342 9.161331 10.75132 0.002702 0.14 0.732792
7.707905 9.24156 10.77522 0.004197 0.16 0.741136
7.833186 9.315465 10.79775 0.00618 0.18 0.748878
7.957108 9.388889 10.82067 0.008916 0.201356 0.756628
8.110717 9.480406 10.85009 0.013738 0.23 0.766383
8.211082 9.540542 10.87 0.017985 0.25 0.772858
8.328215 9.611111 10.894008 0.02431 0.274598 0.78053
8.443105 9.680786 10.91847 0.032232 0.30 0.788191
8.529915 9.733769 10.93762 0.039539 0.32 0.79408
8.655085 9.810744 10.9664 0.052391 0.35 0.802743
8.735663 9.860704 10.98575 0.062285 0.37 0.80844
8.852923 9.934061 11.0152 0.079214 0.40 0.816918
8.92903 9.98214 11.03525 0.091936 0.42 0.822553
9.040583 10.053372 11.06616 0.113227 0.45 0.831024
9.113456 10.100458 11.08746 0.128911 0.47 0.836707
9.214635 10.166667 11.118699 0.153089 0.498239 0.844815
9.29143 10.217655 11.14388 0.17333 0.52 0.851154
9.380601 10.277778 11.174955 0.198872 0.545578 0.858735
9.464812 10.335597 11.20638 0.224964 0.57 0.86613
9.540967 10.388889 11.236811 0.250145 0.592277 0.873034
9.63501 10.456235 11.27746 0.28319 0.62 0.881863
9.702443 10.505738 11.30903 0.308107 0.64 0.888415
9.803071 10.58185 11.36063 0.346963 0.67 0.898554
9.869875 10.634118 11.39836 0.373712 0.69 0.905533
9.969753 10.715403 11.46105 0.414825 0.72 0.916335
10.03617 10.771939 11.50771 0.442722 0.74 0.923751
10.105323 10.833333 11.561344 0.472071 0.760853 0.931646
10.201916 10.924318 11.64672 0.513286 0.79 0.942893
10.268263 10.991105 11.71395 0.54154 0.81 0.95068
10.334807 11.062372 11.78994 0.569668 0.83 0.958435
10.42447 11.166667 11.908864 0.606991 0.856739 0.968563
10.50379 11.268709 12.03363 0.63922 0.88 0.976902
10.57386 11.368284 12.16271 0.666898 0.90 0.983486
10.647043 11.483703 12.32036 0.694863 0.92 0.989289
10.719233 11.611111 12.502989 0.721378 0.938393 0.993718
10.776016 11.722222 12.668428 0.741413 0.951576 0.996241
10.869256 11.928222 12.98719 0.772617 0.97 0.998711
10.933957 12.089833 13.24571 0.79296 0.98 0.9995
11.025138 12.344552 13.66397 0.819725 0.99 0.999907
11.100291 12.577669 14.05505 0.84007 0.995 0.999984
11.176368 12.833333 14.490299 0.859073 0.99781 0.999998
11.239341 13.058332 14.87732 0.873596 0.999 1
11.327018 13.388889 15.45076 0.892028 0.999713 1
11.41125 13.722222 16.033194 0.907834 0.999928 1
11.492504 14.055556 16.618607 0.921391 0.999984 1
11.57158 14.388889 17.206198 0.933078 0.999997 1
11.649017 14.722222 17.795427 0.943174 0.999999 1
11.725194 15.055556 18.385917 0.951895 1 1
11.800381 15.388889 18.977397 0.959419 1 1
11.874778 15.722222 19.569666 0.965894 1 1
11.948534 16.055556 20.162577 0.971449 1 1
12.021764 16.388889 20.756014 0.976199 1 1 ];
%% Matlab GLM fitting
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');
wt = 5 : 0.05 : 17;
[y_fit, dylo, dyhi] = glmval(b, wt, 'probit', stats, 'confidence', 0.95);
mu_1 = -b(1) /b(2);
gradient = central_diff(y_fit, wt);
go = yWT == 1;
nogo = yWT == 0;
fig1 = figure;
plot(xWT(go), yWT(go), 'ko', 'LineWidth', 0.5, 'MarkerFaceColor', 'g', 'MarkerSize', 8)
hold on
plot(xWT(nogo), yWT(nogo), 'ks', 'LineWidth', 0.5, 'MarkerFaceColor', 'r', 'MarkerSize', 8)
p1 = plot(wt, y_fit, '-m');
p2 = plot(wt, y_fit - dylo, '-k');
plot(wt, y_fit + dyhi, '-k')
p3 = plot(wt, gradient/max(gradient), '-', 'Color', 0.6*[1,1,1]);
xlim([5, 17]);
title('Wu and Tian Example')
xlabel('Quantile (q, in units wt)')
ylabel('Probability (p)')
grid on; grid minor
mfw_pos = get(gcf, 'Position'); mfw_pos(3) = mfw_pos(3) * 1.6; set(gcf, 'Position', mfw_pos); % Make Figure Wide (60% wider)
% These are the results from the R code "Gonogo"
% plot(R_output(:,2), R_output(:,5), '--b')
p4 = plot(R_output(:,2), R_output(:,6), '--b');
plot(R_output(:,2), R_output(:,4), '--b')
p5 = plot(R_output(:,1), R_output(:,5), '--r');
plot(R_output(:,3), R_output(:,5), '--r')
legend([p1 p2 p3 p4 p5],{'Logistic Fit','Confidence Bounds from glmval','diff of Logistic Fit','Confidence Interval about p','Confidence Interval about q'}, ...
'Location', 'east')
mdl = fitglm(xWT, yWT, 'Distribution', 'binomial', 'Link', 'probit');
ci = coefCI(mdl, 0.95);
figure(fig1)
plot(wt, mdl.Link.Inverse([ones(size(wt',1),1,"like",wt') wt'] * ci(:,1)), '-c', 'DisplayName', 'from Matlab coefCI')
plot(wt, mdl.Link.Inverse([ones(size(wt',1),1,"like",wt') wt'] * ci(:,2)), '-c', 'HandleVisibility','off')

Réponses (1)

Shubham
Shubham le 25 Juin 2024
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:
  1. Fit the logistic regression model.
  2. Calculate the confidence intervals for the model coefficients.
  3. 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.
  1 commentaire
Craig
Craig le 26 Juin 2024
Hi Shubham,
Thank you very much for the reply! The results that you show perfectly match the results from fitglm and coefCI, plotted in solid cyan and labeled "from Matlab coefCI" in my original plot, but unfortunately do not match the results from Gonogo, plotted in dashed red and labeled "Confidence Interval about q" in my original plot
So, now I know how the fitglm and coefCI method works, but I still don't know how to match the Gonogo results. Does this provide any additional insight for you?

Connectez-vous pour commenter.

Produits


Version

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by