# Logistic regression residual analysis

This question is sort of general and long-winded, but please bear with me.

In my application, I have many datasets, each consisting of ~20,000 datapoints with ~50 features and a single dependent binary variable. I am attempting to model the datasets using regularized logistic regression (R package glmnet)

As part of my analysis, I've created residual plots as follows. For each feature, I sort the datapoints according to the value of that feature, divide the datapoints into 100 buckets, and then compute the average output value and the average prediction value within each bucket. I plot these differences.

Here is an example residual plot:

In the above plot, the feature has a range of [0,1] (with a heavy concentration at 1). As you can see, when the feature value is low, the model appears to be biased towards overestimating the likelihood of a 1-output. For example, in the leftmost bucket, the model overestimates the probability by about 9%.

Armed with this information, I'd like to alter the feature definition in a straightforward manner to roughly correct for this bias. Alterations like replacing

$x \rightarrow \sqrt{x}$

or

$x \rightarrow f_a(x) = \cases{a & if$x<a$\cr x & else}$

How can I do this? I'm looking for a general methodology so that a human could quickly scroll through all ~50 plots and make alterations, and do this for all datasets and repeat often to keep models up-to-date as the data evolves over time.

As a general question, is this even the right approach? Google searches for "logistic regression residual analysis" don't return many results with good practical advice. They seem to be fixated on answering the question, "Is this model a good fit?" and offer various tests like Hosmer-Lemeshow to answer. But I don't care about whether my model is good, I want to know how to make it better!

-

You can't really assess the bias that way in logistic regression. Logisitic regression is only expected to be unbiased on log odds or logit scores, log(p/(1-p)). The proportions will be skewed and therefore look biased. You need to plot the residuals in terms of log odds.

-
How do I combine the log-odd residuals within a bucket? Arithmetic average? This is a little unsettling to me. Intuitively, if a residual analysis shows no bias, then I expect that when the model predicts Pr[y=1]<0.2, then y should equal 1 with probability less than 0.2. But your answer seems to suggest this is not the case. Am I understanding correctly? –  David Jun 26 '12 at 22:30
this is probably better posted as a comment. –  probabilityislogic Jun 26 '12 at 22:47
No David, it doesn't imply anything other than the 0.2 probability, maybe my edits make it more clear. –  John Jun 27 '12 at 3:48
Sorry, I'm still a little confused. My intuitive understanding of an unbiased model is that if the model predicts p=0.2 on every one of a large number of datapoints, then 20% of those datapoints should have y=1. Is this understanding correct? If so, then it seems my plotting methodology should correctly display bias. If not...then I'm not very happy with this concept of "bias"! If an unbiased model reading of 0.2 doesn't tell me anything about the probability that y=1, what good is unbiasedness? –  David Jun 27 '12 at 4:57
Yes, 20% should have y=1. But it's not going to be dead on, it's going to be off by some amount. In the probability space which direction do you think it will be off by and by how much? If it's unbiased it will fall equally in somewhere in the .2:1 or the 0:.2. However, as you can see by the size of those spaces they will tend to be farther away in the bigger area just because they can. In the logit space the distance away should be equal + or -. –  John Jun 27 '12 at 15:13
show 1 more comment

there is unlikely to exist any general software for doing this. most likely because there is no general theory for fixing issues in regression. hence this is more of a "what i would do" type of answer rather than a theoretically grounded procedure.

the plot you produce is basically a visual HL test with 100 bins, but using a single predictor instead of the predicted probability to do the binning. this means your procedure is likely to inherit some of the properties of the HL test.

your procedure sounds reasonable, although you should be aware of "overfitting" your criteria. your criteria is also less useful as a diagnostic because it has become part of the estimation process. also, whenever you do something by intuition, you should write down your decision making process in as much detail as is practical. this is because you may discover the seeds of a general process or theory, which when developed leads to a better procedure (more automatic and optimal with respect to some theory).

i think one way to go is to first reduce the number of plots you need to investigate. one way to do this is to fit each variable as a cubic spline, and then investigate the plots which have non zero non linear estimates. given the number of data points this is also a easy automatic fix for non linearities. this will expand your model from 50 to 200+50k where k is the number of knots. you could think of this as applying a "statistical taylor series expansion" of the "true" transformation.