# How can I compute and plot the 95% confidence bands for a fitted logistic regression model?

I have a set of binary data that I' d like to fit a logistic regression model to and show the 95% confidence band around the fitted function. I can produce the regression model but how do I compute and show the confidence bands in plot3 ?

If showing the confidence bands is not possible can the single prediction bands be shown ? Many thanks.

damagecodedata = {{53, 1}, {57, 1}, {58, 1}, {63, 1}, {66, 0}, {67, 0},
{67, 0}, {67, 0}, {68, 0}, {69, 0}, {70, 1}, {70, 0},
{70, 1}, {70, 0}, {72, 0}, {73, 0}, {75, 1}, {75, 0},
{76, 0}, {76, 0}, {78, 0}, {79, 0}, {81, 0}};

plot1 = ListPlot[damagecodedata, PlotRange -> {{50, 90}, {-.1, 1.1}},
Frame -> True,
FrameStyle -> Directive[FontSize -> 13],
FrameLabel -> {Style["Temperature (deg F)", 14], Style["Coded Damage", 14]},
FrameTicks -> {Range[50, 90, 5], {0, 1}, None, None},
PlotStyle -> {PointSize[Medium], Red}, ImageSize -> 400
]

model = LogitModelFit[damagecodedata, x, x]
Normal[model]
model["ParameterTable"]

plot2 = Plot[model[x], {x, 50, 90}];
plot3 = Show[plot1, plot2] • The documentation for LogitModelfit has information on confidence intervals. Does model["ParameterConfidenceIntervals"] not give what you want? Jun 7, 2013 at 14:33
• Doesn't ParameterConfidenceIntervals give us the intervals for the 2 parameters only {{0.581069, 29.5047}, {-0.444302, -0.0200235}} ? I'd like to show the confidence bands for the entire model plotted in plot3 on either side of the blue fitted curve. Jun 7, 2013 at 14:38
• I'm horrible with stats, but I'm assuming if they're not included in the model parameters, it's likely because they're undefined for this model. Jun 7, 2013 at 14:47
• In the documentation for NonlinearModelFit there is an example under Applications (Obtain and visualize 90% confidence bands for the fit). This is what I'd like to produce but using LogitModelFit. Jun 7, 2013 at 14:49
• Corey, good point, I'm not a statistician either. If what you say is true, then can anyone offer how I could use the ParameterConfidenceIntervals info. to produce a visual performance indicator for the plot3 curve fit ? Jun 7, 2013 at 14:53

It does not appear that LogitModelFit and GeneralizedLinearModelFit have the options to produce confidence bands for predictions but those are easily constructed.

(* Get parameter estimates and the estimate of the covariance matrix *)
parms = model["BestFitParameters"];
cov = model["CovarianceMatrix"];

(* Generate a list of temperatures to make predictions *)
tMin = Min[damagecodedata[[All, 1]]];
tMax = Max[damagecodedata[[All, 1]]];
t = Table[tMin + (tMax - tMin) i/100, {i, 0, 100}];

(* Get predictions and upper and lower confidence limits for a + b*t *)
pred = Table[1 - 1/(1 + Exp[parms.{1, t[[i]]}]), {i, 101}];
lower = Table[1 - 1/(1 + Exp[parms.{1, t[[i]]} - 1.96 Sqrt[{1, t[[i]]}.cov.{1, t[[i]]}]]), {i, 101}];
upper = Table[1 - 1/(1 + Exp[parms.{1, t[[i]]} + 1.96 Sqrt[{1, t[[i]]}.cov.{1, t[[i]]}]]), {i, 101}];

(* Plot results *)
ListPlot[{Transpose[{t, pred}], Transpose[{t, lower}], Transpose[{t, upper}]},
Joined -> True, PlotStyle -> {Black, Gray, Gray}, Frame -> True,
FrameLabel -> {"Temperature (deg F)", "Probability"}] What you are getting is a set of 95% confidence intervals for the parameter $1-1/(1+e^{a+b t})$ (the probability of a success - damage or no damage?) for a range of values of temperature (t).

Update

One can do this much shorter by using Plot rather than ListPlot:

(* Get parameter estimates and the estimate of the covariance matrix *)
parms = model["BestFitParameters"];
cov = model["CovarianceMatrix"];

(* Plot results *)
Plot[{1 - 1/(1 + Exp[parms.{1, t}]),
1 - 1/(1 + Exp[parms.{1, t} - 1.96 Sqrt[{1, t}.cov.{1, t}]]),
1 - 1/(1 + Exp[parms.{1, t} + 1.96 Sqrt[{1, t}.cov.{1, t}]])},
{t, Min[damagecodedata[[All, 1]]], Max[damagecodedata[[All, 1]]]},
PlotStyle -> {Black, Gray, Gray},
Frame -> True, FrameLabel -> {"Temperature (deg F)", "Probability"}]


But then you wouldn't get a list of values if you needed to put those into a table in a report.

• Can this be extended to more than one independent variable (x1, x2, ...)? Even more so with interpendency ( bx*x1*x2)? May 14, 2021 at 15:31
• @Eisbär. Certainly. For a similar figure you'd just fix the other predictor variables at some value. But a 3D "volume" figure could also be constructed if another predictor is allowed to vary.
– JimB
May 14, 2021 at 15:36
• Precisely answered! Would you answer an according post I would like to create? May 14, 2021 at 15:39
• @Eisbär I'll definitely take a look at it. And there are many others much better at this than me.
– JimB
May 14, 2021 at 15:41
• mathematica.stackexchange.com/questions/246038/… May 14, 2021 at 16:07

I'm no expert in statistics, and what I've written below is based on some reading from this resource.

The Logistic Regression Model

The more familiar(?) linear regression model is typically used to analyze the relationship between a quantitative response variable and an explanatory variable. In the logistic regression model, the response variable has only two possibilities. Take for example, the flipping of a coin. There are only two possibilities (heads or tails) and the chance of the flip result being heads can be modeled using a binomial distribution. However, if one were to ask if the height of the flip influences the outcome, we now have an explanatory variable (the height) and we use the logistic regression model to explore the effect of height on the outcome of the flip.

Comparison with linear regression

In a linear regression, we model the mean (u) of the response variable (y) as a linear function of the explanatory variable: u = b0 + b1 * x. In a logistic model, rather than working with means, we work with odds, p, (e.g. the odds of a coin flip resulting in heads). One could substitute p for u in the linear expression above; however one would obtain a probability greater than 1 for large values of x if b1 is not equal to 0. To avert this problem, the logistic model uses log(p) rather than p, giving us a model:

log(p/(1-p))= b0 + b1 x

where p is the binomial proportion, x is the explanatory variable and the model parameters are b0 and b1.

The special case of a binary explanatory variable

The reference above gives an adequate description of the case of a binary explanatory variable, which is useful in that (a) it can be calculated by hand and (b) it helps to understand the significance of the model parameters. Rather than rehash the reference, allow me to say that the slope (b1) is the odds ratio which is the increase in odds of an acceptable outcome by increasing the explanatory variable by 1 unit.

Here we come (I think) to the crux of the question, which is how one plots the confidence interval of the model. As I understand the logistic model, plotting a confidence interval is not providing meaningful information since the confidence interval is the range in which there is (in this case) a 95% chance of finding the actual slope.

My attempt at analyzing the data

For sake of simplicity (giving an odds ratio > 1) I am going to invert the OP coding

damagecodedata = {{53, 1}, {57, 1}, {58, 1}, {63, 1}, {66, 0}, {67, 0},
{67, 0}, {67, 0}, {68, 0}, {69, 0}, {70, 1}, {70, 0},
{70, 1}, {70, 0}, {72, 0}, {73, 0}, {75, 1}, {75, 0},
{76, 0}, {76, 0}, {78, 0}, {79, 0}, {81, 0}};
damagecodedata[[All,2]] = damagecodedata[[All,2]]/.{0->1,1->0};
model = LogitModelFit[damagecodedata, x, x]
Normal[model]
model["ParameterTable"]


The model table gives a slope of 0.232 with an error of 0.108, both of which are in log form. The slope and standard error are extracted:

model["ParameterTableEntries"][[2, {1, 2}]]


and we obtain the odds ratio with confidence interval:

E^{a, a + b, a - b} /. {a -> model["ParameterTableEntries"][[2, 1]],
b -> 1.96*model["ParameterTableEntries"][[2, 2]]}


yielding an odds ratio of 1.26 with a range of 1.02 to 1.55. Meaning that the chance of a (1) outcome (remember, I changed the definition of the outcome from the OP) increases by 1.26 times with every 1 unit change in temperature.

Yeah, but what about the graph!

The OP refers to page 14 of this link as an example of the graph that is desired. There is a fundamental difference in the data in the reference and the data presented here, namely the reference is using probabilities that either aren't known or haven't been presented in this question. If the data are available and it is possible to use probabilities instead of the boolean response, then one might be able to generate the requested graph. Presently, I am not sure that the graph would provide any meaningful information.

Update

We can use the "ParameterConfidenceIntervals" property of LogicModelFit to generate confidence bands, although to be honest I do not know what their value is for a logistic model in general and as you will see below, they are not at all useful to the current problem. Nonetheless, I'm presenting how the confidence bands can be drawn using Mathematica

To make the plotting easy, I make a list of parameters and define the logistic function.

res = Join[{model["BestFitParameters"]}, model["ParameterConfidenceIntervals"] // Transpose]
f[{a_, b_}] := 1/(1 + E^(-a - b x));


Then we can Plot:

Plot[Evaluate[f[#] & /@ res], {x, 53, 81}, Filling -> {{3 -> {1}}, {2 -> {1}}},
FillingStyle -> {LightBlue, LightGreen}, PlotRange -> All, Frame -> True,
FrameStyle -> Directive[FontSize -> 13], FrameLabel ->
{Style["Temperature (deg F)", 14], Style["Coded Damage", 14]},
FrameTicks -> {Range[50, 90, 5], {0, 1}, None, None},
Epilog -> {PointSize[Medium], Red, Point@damagecodedata}]


Giving us the following: As can be seen, the 95% confidence intervals as presented from the model fit tell us that the model isn't very good for the data set.

• thank you for your thoughtful response. I will reply again after I've digested what you have written. I'm realizing that my question has as much to do with probability theory as Mathematica usage and may not have phrased the question as precisely as possible. Upon reflection I'm really trying to determine how trustworthy the blue curve in plot3 is. Everything else equal, if we had 230 points instead of 23 we would expect the level of trustworthiness to be higher and with only 5 points we would have less trust. Maybe the question is "how trustworthy is the blue curve "? Jun 7, 2013 at 17:48
• One thing to keep in mind (and as I'm not a statistician, I'm probably using the wrong terminology here) is that there are two types of uncertainty to discuss: one is the uncertainty in the slope (how confident are you in the increase in odds with increase in explanatory variable) and how confident you are in the probability of an outcome at a given value of the explanatory variable. As I understand it, you are presently more interested in this latter value, which will require additional data points acquired at each temperature. Jun 7, 2013 at 18:07
• Yes, I'm more interested in the confidence of the damage code as a function of temp. Let's drill a little deeper for a partial answer. Two important points on the fitted curve are (52F,0.95) and (77F, 0.05). The model is saying for temps less than 52F there is a better than 95% prob of damage. And for temps greater than 77F there is a better than 95% prob of no damage. In between these 2 boundaries we have higher uncertainty (per the nominal model). I wonder if the problem would be easier if Mathematica could be focused on the "tails" of the function only ? The residuals also play a role here. Jun 7, 2013 at 18:54
• I suspect what we are running up against is the validity of the model. Your data look very similar to the Challenger Distaster data found here: stat.wisc.edu/~mchung/teaching/MIA/reading/… and the point made on page 10 of that document, namely that logistic transformations require most of the data not being 1's and 0's is adding to the difficulty of your data analysis. Jun 7, 2013 at 19:38
• This is exactly my data (nice detective work). I didn't want answers to simply point me to work already done. Plus my inquiry is independent of the data source and others have confirmed my parameter values removing the need to question the parameter values. So the question still remains if for example the temperature is say, 52F and the model says there is a 95% chance of damage what is the tolerance, confidence or degree of trustworthiness in the model prediction ? Is it +5%/-60% or +1%/-2% . . . ? There must be an answer to this. Every model has uncertainty associated with it. Jun 7, 2013 at 19:56

I'm wondering if you ever got this worked out, because I'm trying to do the same thing.

Here's how far I got. You can retrieve confidence intervals for the fitted curve if you use NonLinearModelFit rather than LogitModelFit. The syntax for a logistic function in NonLinearModelFit is:

model2 = NonlinearModelFit[damagecodedata, (1/(1 + E^(-a + b*x))),
{{a, -15}, {b, 0.1}}, x]


Note I used initial guesses for the parameters that were reasonably close to what you got with LogitModelFit. The fitted model is:

1/(1 + E^(-19.2615 + 0.300954 x))


The parameters are a little different from what you got, because NonlinearModelFit assumes the data are normally distributed whereas LogitModelFit assumes a binomial distribution (I think the latter is more reasonable).

You can obtain the confidence bands like so:

model2["MeanPredictionBands", ConfidenceLevel -> .90]


And when you plot them you get: The problem, as I see it, is the confidence bands go above 1 and below 0, even though the function cannot.

Any help with this would be appreciated.

• Not exactly an answer to the question at hand -- more of a (semi) workaround. You can comment directly below a post as the answer section is, well, strictly for answers. Mar 28, 2016 at 19:55
• @Sektor That's a bit harsh, isn't it? I quite often see partial workarounds as answers. Besides, this answer/comment, with its graphics and code, wouldn't fit in a comment anyway. Mar 28, 2016 at 20:34
• BTW I ended up deciding to calculate the confidence bands for my particular logistic regression problem by bootstrapping. I sampled my original array of data 1000 times with replacement, using the command sample = RandomChoice[data, Length[data]] Mar 29, 2016 at 4:54
• [cont'd] Then I calculated a logistic fit for each sample with LogitModelFit. Then for an array of 50 x-axis values, I calculated the fitted function values for each of the 1000 fitted logistic functions. I sorted them low to high and took the 100th- and 900th-highest values to be my desired confidence interval. As one would hope, the confidence band never goes above y=1 or below y=0. Mar 29, 2016 at 5:06
• @SjoerdC.deVries Well, I certainly didn't want it to sound harsh. Decided to write a comment myself instead of macro'ing the guy. Plus, I said myself that it's a workaround, which by definition is a solution to the problem. Nevetheless it was a suggestion to improve on the answer Mar 29, 2016 at 8:02

Imagine your fitted model is called nlm (as named in the documentation center examples).

Then, throguh the command:

nlm["MeanPredictionConfidenceIntervals", ConfidenceLevel -> 0.95]


you can obtain the confidence interval band for your model. In this case I consider the bands at 95% Condifence Level, you can instead use 0.682689492 to obtain the bands at one sigma.

• I tried your suggestion and it simply returned my input to me. Moreover, if you execute model["Properties"] MeanPredictionConfidenceIntervals doesn't show up on the list. So I don't think MeanPredictionConfidenceIntervals is a valid property for LogitModelFit. If you think otherwise can you post your code (inputs and outputs) so I can confirm your results ? Thanks. Oct 9, 2013 at 16:47