# Binary (binomial) logistic regression with multiple independent variables: calculating the confidence range or prediction

Adding on top to this post here, I am looking for the confidence range (intreval) of the prediction of a multi variant binary logistic regression.

Can someone extend the solution in the post referenced above or provide an alternative?

(*Test data (not full set)*)
testDate = {{449., 30., 0.}, {449., 30., 0.}, {452., 30.,
0.}, {452., 30., 0.}, {562., 30., 0.}, {566., 30.,
1.}, {566., 30., 1.}, {656., 30., 1.}, {656., 30.,
1.}, {658., 30., 1.}, {658., 30., 1.}, {452., 34.,
0.}, {452., 34., 0.}, {452., 34., 0.}, {562., 34.,
0.}, {562., 34., 0.}, {562., 34., 0.}, {658., 34.,
0.}, {658., 34., 1.}, {658., 34., 1.}, {449., 50.,
1.}, {449., 50., 1.}, {449., 50., 1.}, {449., 50.,
1.}, {452., 50., 1.}, {562., 50., 1.}, {562., 50.,
1.}, {566., 50., 1.}, {566., 50., 1.}, {566., 50.,
1.}, {566., 50., 1.}, {656., 50., 1.}, {656., 50.,
1.}, {656., 50., 1.}, {658., 50., 1.}, {658., 50., 1.}};
(*Model to fit*)
logit = LogitModelFit[testDate, {P, p, P*p}, {P, p}];
Normal[logit]


Here is the basic approach to get confidence intervals for the probability associated with some predictor values.

testData = {{449, 30, 0}, {449, 30, 0}, {452, 30, 0}, {452, 30, 0}, {562, 30, 0}, {566, 30, 1},
{566, 30, 1}, {656, 30, 1}, {656, 30, 1}, {658, 30, 1}, {658, 30, 1}, {452, 34, 0}, {452, 34, 0},
{452, 34, 0}, {562, 34, 0}, {562, 34, 0}, {562, 34, 0}, {658, 34, 0}, {658, 34, 1}, {658, 34, 1},
{449, 50, 1}, {449, 50, 1}, {449, 50, 1}, {449, 50, 1}, {452, 50, 1}, {562, 50, 1}, {562, 50, 1},
{566, 50, 1}, {566, 50, 1}, {566, 50, 1}, {566, 50, 1}, {656, 50, 1}, {656, 50, 1}, {656, 50, 1},
{658, 50, 1}, {658, 50, 1}};
(*Model to fit*)
logit = LogitModelFit[testData, {P, p, P*p}, {P, p}];

(* Get the necessary information from the logistic fit *)
parms = logit["BestFitParameters"]
cov = logit["CovarianceMatrix"]

(* Generate data points for which you want predictions and confidence intervals *)
(* Here we'll start out with just making predictions for the values of the predictors
in the original dataset.  I'm assuming that a "1" isn't included already for the
intercept. *)
predictors = {1, #[[1]], #[[2]], #[[1]] #[[2]]} & /@ testData

(* Make predictions *)
predictions = 1 - 1/(1 + Exp[predictors . parms])
t = 1.96;  (* t=InverseCDF[StudentTDistribution[Length[testData]-Length[parms]],0.975]
if you're conservative or want to match what SAS does *)
lower = 1 - 1/(1 + Exp[predictors . parms - t Diagonal[predictors . cov . Transpose[predictors]]^0.5])
upper = 1 - 1/(1 + Exp[predictors . parms + t Diagonal[predictors . cov . Transpose[predictors]]^0.5])


Using Manipulate can be helpful to see the effect of one or more variables on a range of values for a particular predictor variable. (But "P" and "p" aren't visually very separate labels for this example.)

Manipulate[pMin = Min[testData[[All, 2]]];
pMax = Max[testData[[All, 2]]];
predictors = {1, P, p, P*p};
predictions = 1 - 1/(1 + Exp[predictors . parms]);
t = 1.96;
lower = 1 - 1/(1 + Exp[predictors . parms - t (predictors . cov . Transpose[predictors])^0.5]);
upper = 1 - 1/(1 + Exp[predictors . parms + t (predictors . cov . Transpose[predictors])^0.5]);
Plot[{predictions, lower, upper}, {p, pMin, pMax},  PlotRange -> {Automatic, {0, 1}},
PlotStyle -> {Black, {Black, Dotted}, {Black, Dotted}},
Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"p", "Probability"}],
{{P, Mean[testData[[All, 1]]]}, Min[testData[[All, 1]]],
Max[testData[[All, 1]]], Appearance -> "Labeled"}]


When there are two predictors you can also show a 3D plot:

Plot3D[{predictions, lower, upper}, {P, Min[testData[[All, 1]]], Max[testData[[All, 1]]]},
{p, Min[testData[[All, 2]]], Max[testData[[All, 2]]]}, BoxRatios -> {1, 1, 1},
PlotStyle -> {{Blue, Opacity[1]}, {Cyan, Opacity[0.5]}, {Cyan, Opacity[0.5]}},
AxesLabel -> (Style[#, Bold, 18] &) /@ {"P", "p", "Probability"}, ImageSize -> Large]


• For me the Manipulate only works if I take out the Transpose in lower and upper. Apart from that it works greatly. May 15, 2021 at 17:10
• Thanks for the Accept. I wish I could explain the issue with Transpose but with a fresh kernel on Windows 10 Mathematica 12.2, it works fine as is for me.
– JimB
May 15, 2021 at 18:41
• predictors in Manipulate is a list. This cannot be transposed. (e.g. Transpose[{1,2,3}] ). May 15, 2021 at 19:06
• Very odd. Do you have a different version of Mathematica? It works for me with and without Transpose.
– JimB
May 15, 2021 at 19:11
• I am at 12.1.1.0 version. Anyhow, it seems to not care anymore if there is a list to be transposed in 12.2. Good time to update. Btw. I referenced the code here with Minitab result. It matches! Nice programming and thank you for the work put into this answer. Can you give a reference for the method? May 15, 2021 at 19:16