# 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? – bobthechemist Jun 7 '13 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. – Steve Jun 7 '13 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. – Corey Kelly Jun 7 '13 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. – Steve Jun 7 '13 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 ? – Steve Jun 7 '13 at 14:53

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 "? – Steve Jun 7 '13 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. – bobthechemist Jun 7 '13 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. – Steve Jun 7 '13 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. – bobthechemist Jun 7 '13 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. – Steve Jun 7 '13 at 19:56

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. – Steve Oct 9 '13 at 16:47