Measuring accuracy of a logistic regression-based model

Question

I have a trained logistic regression model that I am applying to a testing data set. The dependent variable is binary (boolean). For each sample in the testing data set, I apply the logistic regression model to generates a % probability that the dependent variable will be true. Then I record whether the acutal value was true or false. I'm trying to calculate an $R^2$ or Adjusted $R^2$ figure as in a linear regression model.

This gives me a record for each sample in the testing set like:

prob_value_is_true         acutal_value
   .34                          0
   .45                          1
   .11                          0
   .84                          0
    ....                        ....          

I am wondering how to test the accuracy of the model. My first attempt was to use a contingency table and say "if prob_value_is_true > 0.80, guess that the actual value is true" and then measure the ratio of correct to incorrect classifications. But I don't like that, because it feels more like I'm just evaluating the 0.80 as a boundary, not the accuracy of the model as a whole and at all prob_value_is_true values.

Then I tried to just look at each prob_value_is_true discrete value, as an example, looking at all samples where prob_value_is_true=0.34 and measuring the % of those samples where the acutal value is true (in this case, perfect accuracy would be if the % of samples that was true = 34%). I might create a model accuracy score by summing the difference at each discrete value of prob_value_is_true. But sample sizes are a huge concern here, especially for the extremes (nearing 0% or 100%), such that the averages of the acutal values are not accurate, so using them to measure the model accuracy doesn't seem right.

I even tried creating huge ranges to ensure sufficient sample sizes (0-.25, .25-.50, .50-.75, .75-1.0), but how to measure "goodness" of that % of actual value stumps me. Say all samples where prob_value_is_true is between 0.25 and 0.50 have an average acutal_value of 0.45. Is that good since its in the range? Bad since its not near 37.5% (the center of the range)?

So I'm stuck at what seems like should be an easy question, and hoping someone can point me to a resource or method to calculate an accuracy stastic for a logistic regression model.

Answer 1

A measure that is often used to validate logistic regression, is the AUC of the ROC curve (plot of sensitivity against 1-specificity - just google for the terms if needed). This, in essence, evaluates the whole range of threshold values.

On the downside: evaluating the whole range of threshold values may be not what you're after, since this (typically) includes thresholds that result in very large numbers of false negatives or false positives. There are versions of the AUC that account for this (partial AUC), so if that is an issue for you, you can look into that.

Answer 2

You are correct to worry about proportion classified correct as mainly reflecting the effect of an arbitrary boundary. I'd recommend two measures. One is the $c$-index or ROC area as others have described. This has an interpretation that is simpler than thinking about an ROC curve, and is a measure of pure predictive discrimination. Secondly, estimate a continuous calibration curve without any binning of data. If predictions are being assessed on an independent dataset, you can use lowess with outlier detection turned off to estimate the relationship between predicted and actual Prob[Y=1]. The val.prob function in the R rms package will do both of these things. Other functions in rms will do the same for internal validation, using resampling to remove the effects of overfitting.

Answer 3

If your data are grouped by $x$ values, you can compute the model predicted value and it's associated confidence interval, and see if the observed percentage falls within that range. For example, if you had 10 observations at $x=10$, 10 obs at $x=20$, 10 obs at $x=30$, etc., then mean(y[x==10]==1), mean(y[x==20]==1), etc., would yield percentages that can be compared to predictions. Bear in mind, that even if the model is perfect, some observed percentages will bounce outside of the 95% CI, just like in OLS regression. If your data are not grouped, you can form your own groups by binning the data according to ranges of the $x$ variable, as you suggest. This isn't fully valid, as it will depend on the choice of bins, can be useful as a way of exploring your model.

In general, the task you have given yourself here is difficult. That's because, with logistic regression, you are dealing with two different kinds of things. The model's predictions are a latent variable, whereas your observed response variable (while presumably generated by a latent variable) is not. Of course, people will often want to know what the predicted response is, and that's totally reasonable; this is just one of those cases where life isn't fair.

If you do want to predict the outcome, you need to decide what you want to maximize. If you have just 1 case, and you want your prediction to be most likely to be right, you should predict $y=1$, if $\hat y\ge .5$. (This is all pretty intuitive.) On the other hand, if you want to maximize overall accuracy over your total sample (or any other group), you should predict $y=1$, if $\hat y \ge p(y=1)$. For example, let's say that in your sample, 30% of all cases are 1's, then if $\hat y = .31$, you should predict that $y$ will be $1$, even though it's $<.5$. This is counter-intuitive, and a lot of people stumble here, but this algorithm will maximize your accuracy.

A more comprehensive way to think about how much information is in your model, is to integrate over how accurate you would be given every possible threshold $(0, 1)$. This is the area under the curve (AUC) of the model's receiver operating characteristic (ROC), discussed by @Nick Sabbe. Remember that there is no $R^2$ for logistic regression. There are so called 'pseudo $R^2$'s, but the AUC (or the concordance, $c$, a synonym) is probably the best way to think about this issue.

Answer 4

You may want to have look at my package softclassval (at softclassval.r-forge.r-project.org you find also two oral presentations I gave about the ideas behind the package).

I wrote it for a slightly different problem, namely if the reference (e.g. pathologist) "refuses" to give a clear class. However, you can use it with "normal" hard classes and it avoids the definition of a threshold for hardening the originally continuous prediction - so you don't evaluate the 0.8.

However, I recommend to use it alongside with, say, a ROC or specificity-sensitivity-diagram: the results will often look quite bad as "my" methods will penalize already slight deviations (e.g. 0.9 instead of 1 gives 0.1 difference for my measures, but all thresholds below 0.9 will ignore this). Actually I think is rather an advantage: the lack of this sensitivity agaist small deviations is one of the major points of criticism with those "hardened" measures like accuracy, sensitivity, recall, etc.

In addition, by comparing mean absolute error (MAE) and root mean squared error RMSE you can find out whether you have many small deviations or fewer grossly misjudged samples.

Answer 5

I wonder why you aren't using the bernoulli log-likelihood function. Basically, for every $0$ actual value, you score $-\log (1-\hat {p}) $. This measures how close to predicting $0$ your model is. Similarly, for every $1$ actual value you score $-\log (\hat {p}) $. This measures how close to predicting $1$ your model is.

This does not suffer from arbitrary thresholds. The smaller the measure the better.

Answer 6

I think you could establish a threshold (say 0.5), so when your probability is equals to or greater than that threshold your predicted class would be 1, and 0 otherwise. Then, you could obtain a measure of your accuracy in this way:

confusion_matrix <- ftable(actual_value, predicted_value)
accuracy <- sum(diag(confusion_matrix))/number of events*100

Given that your probability is the probability of given your data (x) and using your model your class value (y) is equal to 1, I do not understand why you always obtain probability values lower than 0.5. Which is the frequency of your actual classes (actual_value)?

Answer 7

Here's my quick suggestion: Since your dependent variable is binary, you can assume it follows a Bernoulli distribution, with probability given by logistic regression $Pr_{i} = invlogit(a + bx_{i})$.

Now, set one simulation as follow: $ y.rep[i] \sim Bernoulli (p[i])$

Then, run this simulation, say, 100 times. You will have a matrix with n rows (n is the number of subjects) and k columns (in this case, k=100, the number of simulations). In r code:

for (j  in 1:100)
  mat.y.rep[,j] <- Bernoulli ( p) # p is a vector with a probability for each subject

Now you compute the difference between the predicted in each simulation and observed. After computing this difference, just compute the mean number of true-positive and false-positive for each row (each subject) and plot the histogram. Or compute both for each column (simulation) e plot the histogram (I prefer this).

Hope it helps...

Answer 8

There are many ways to estimate the accuracy of such predictions and the optimal choice really depends on what will the estimation implemented for.

For example, if you plan to select a few high score hits for an expensive follow-up study you may want to maximize the precision at high scores. On the other hand, if the follow-up study is cheap you may want to maximize the recall (sensitivity) at lower scores. The ROC AUC may be suitable if you are comparing different method, etc.

On the practical side, R's ROCR package contains 2 useful functions

pred.obj <- prediction(predictions, labels,...)
performance(pred.obj, measure, ...)

Together, these functions can calculate a wide range of accuracy measures, including global scalar values (such as "auc") and score-dependent vectors for plotting Recall-precision and ROC curves ("prec", "rec", "tpr" and "fpr", etc.)

Answer 9

You need to define what you mean by "accuracy". What you would like to know, please pardon me for putting words in your mouth, is how well your model fits the training data, and more importantly, how well this model "generalizes" to samples not in your training data. Although ROC curves can be useful in analyzing the tradeoff between precision and recall for various values of the threshold, I suggest adding mean-squared-error, or the Brier score to your toolbox. It's easy to compute, and you can immediately get a feel for whether feature changes affect the fit of the model, when applied to training data. Since overfit is possible in this case, your job isn't done here. To evaluate generalization performance, or how well you do on data you haven't seen, it isn't enough to look at your performance on the training samples. Of course your model is good at those, because they the values you used to determine the coefficients to your logistic. You need to set aside some samples for test data. Your MSE performance on this set should set your generalization expectations according to the Hoeffding inequality. Your maximum generalization error will depend on the number of features in your model as well as the number of samples used to compute the test statistic. Be mindful that you'll need to steal some of your training samples for test samples. I recommend 10-fold cross-validation, where you shuffle, choose 90% for training, 10% for testing, and then measure, repeat, and then average all the measurements.