4
$\begingroup$

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]
$\endgroup$

1 Answer 1

5
$\begingroup$

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"}]

Manipulate example of p vs probability for various values of P

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]

Prediction surface and 95% confidence limits

$\endgroup$
5
  • $\begingroup$ For me the Manipulate only works if I take out the Transpose in lower and upper. Apart from that it works greatly. $\endgroup$
    – Eisbär
    May 15, 2021 at 17:10
  • $\begingroup$ 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. $\endgroup$
    – JimB
    May 15, 2021 at 18:41
  • $\begingroup$ predictors in Manipulate is a list. This cannot be transposed. (e.g. Transpose[{1,2,3}] ). $\endgroup$
    – Eisbär
    May 15, 2021 at 19:06
  • $\begingroup$ Very odd. Do you have a different version of Mathematica? It works for me with and without Transpose. $\endgroup$
    – JimB
    May 15, 2021 at 19:11
  • $\begingroup$ 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? $\endgroup$
    – Eisbär
    May 15, 2021 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.