5
$\begingroup$

Is the probability calculated by a logistic regression model (the one that is logit transformed) the fit of cumulative distribution function of successes of original data (ordered by the X variable)?

EDIT: In other words - how to plot the probability distribution of the original data that you get when you fit a logistic regression model?

The motivation for the question was Jeff Leak's example of regression on the Raven's score in a game and whether they won or not (from Coursera's Data Analysis course). Admittedly, the problem is artificial (see @FrankHarrell's comment below). Here is his data with a mix of his and my code:

download.file("http://dl.dropbox.com/u/7710864/data/ravensData.rda", 
              destfile="ravensData.rda", method="internal")
load("ravensData.rda")

plot(ravenWinNum~ravenScore, data=ravensData)

enter image description here

It doesn't seem like good material for logistic regression, but let's try anyway:

logRegRavens <- glm(ravenWinNum ~ ravenScore, data=ravensData, family=binomial)
summary(logRegRavens)
# the beta is not significant

# sort table by ravenScore (X)
rav2 = ravensData[order(ravensData$ravenScore), ]

# plot CDF
plot(sort(ravensData$ravenScore), cumsum(rav2$ravenWinNum)/sum(rav2$ravenWinNum), 
     pch=19, col="blue", xlab="Score", ylab="Prob Ravens Win", ylim=c(0,1), 
     xlim=c(-10,50))
# overplot fitted values (Jeff's)
points(ravensData$ravenScore, logRegRavens$fitted, pch=19, col="red")
# overplot regression curve
curve(1/(1+exp(-(logRegRavens$coef[1]+logRegRavens$coef[2]*x))), -10, 50, add=T)

If I understand logistic regression correctly, R does a pretty bad job at finding the right coefficients in this case.

enter image description here

  • blue = original data to be fitted, I believe (CDF)
  • red = prediction from the model (fitted data = projection of original data onto regression curve)

SOLVED
- lowess seems to be a good non-parametric estimator of the original data = what is being fitted (thanks @gung). Seeing it allows us to choose the right model, which in this case would be adding squared term to the previous model (@gung)
- Of course, the problem is pretty artificial and modelling it rather pointless in general (@FrankHarrell)
- in regular logistic regression it's not CDF, but point probabilities - first pointed out by @FrankHarrell; also my embarrassing inability to calculate CDF pointed out by @gung.

$\endgroup$
2
$\begingroup$

There are several issues here. The biggest is that the thinking behind the code for generating the blue points in the second plot is confused. I can't make any sense of what the logic is supposed to be. The $x$ values are the same as your data, but the $y$ values are the number of wins up to that point on the $x$ axis divided by the total number of wins. This really makes no sense. You aren't even dividing by the total number of games, which means you don't have the cumulative probability at that point (although that still wouldn't make sense). This also makes you wonder why the last blue point isn't at $1.0$. Actually it is, but the last point is at $x=55$, which is excluded by your argument xlim=c(-10,50).

At any rate, if you want to get a non-parametric estimate of what the underlying function might be, you could plot a lowess line:

plot(ravenWinNum~ravenScore, data=ravensData)
lines(with(ravensData, lowess(ravenWinNum~ravenScore)))

enter image description here

What is clear from this figure, is that there is a curvilinear relationship between the score and the probability of winning (at least in these data, we can debate whether that makes any theoretical sense). If you thought about it, it was clear from the plot of the raw data in that there are increasing numbers of wins at both the low end and high end of $x$; the losses are in the center. A very simple way to fit a curvilinear relationship is to add a squared term:

lr.Ravens2 <- glm(ravenWinNum~ravenScore+I(ravenScore^2), data=ravensData, 
                  family=binomial)

We can test these models as a whole by assessing their deviance against the null deviance and comparing the result to a chi-squared distribution with the number of degrees of freedom consumed. We do this because we want to think of both the linear term and the squared term together as a unit. This method of testing the model as a whole is analogous to the global $F$-test that comes with a multiple regression model.

# here is the test of the model with the squared term as a whole
1-pchisq(24.435-16.875, df=2)
[1] 0.02282269
# here is the test of your model
1-pchisq(24.435-20.895, df=1)
[1] 0.05990546
# this is the test of the improvement due to adding the squared term
1-pchisq(20.895-16.875, df=1)
[1] 0.04496371

I wouldn't make too much of the arbitrary threshold of significance at $.05$, but we see that the initial model is not significant by conventional standards, or at least is slightly less significant. (Given that you have so few data, I would think a higher $\alpha$ is fine.) The model that includes the squared term as well is more significant, even though it uses up another of your precious few degrees of freedom. The AIC is lower as well: the first model has $AIC=24.895$, whereas the second model has $AIC=22.875$. The last of the tests listed shows that the addition of the squared term yields a significant improvement in fit.

The implication of the fact that there is a curvilinear relationship between $x$ and $y$ is that the function that maps the conditional probability that $y_i=1$ given $x_i$ cannot be equated to a cumulative probability distribution that takes $x$ as its quantiles.

For completeness, here is the same plot with the model's predicted probabilities overlaid:

enter image description here

We see that the predicted probabilities mirror the lowess line reasonably well. Of course, they don't match perfectly, but there is no reason to believe that the lowess line is exactly correct; it just provides a non-parametric approximation of whatever the underlying function is.


Having addressed the specific example in your question, it is worth addressing the assumption underlying the question more broadly. What is the relationship between logistic regression and a cumulative distribution function (CDF)?

There can be a connection between binomial regression and a CDF. For example, when the probit link is used, the conditional probability of 'success' ($\pi_i(Y_i=1|X=x_i)$) is transformed via the inverse of the normal CDF. (For more about that, it may help to read my answer here: difference between logit and probit models.) But there are a couple of points to be made here. First, as @FrankHarrell and @Glen_b note, binomial (e.g., logistic) regression is about modeling conditional probabilities at each point.

Second, although the fitted function can look sort of like a CDF when the conditions are just right, it will never actually be a function (scaled, shifted, etc.) of the CDF of your $X$ data. Your own data serve as an example of the function not even looking like a CDF because you have a polynomial relationship. But it is true even when the function does look like a CDF. The easiest way to understand this is that your last (highest $x$ valued) point will always be 100% of the way through your data, but the fitted function will never reach $\hat p_i=1.0$, it can only approach that value as $X$ (the variable, not your data) approaches infinity. As a result, you should not think of a function of the CDF of your $x$ data as the values the logistic regression is supposed to approximate.

$\endgroup$
  • $\begingroup$ You hit the nail on the head! First - lowess. I started thinking moving average (which technically is lowess of 0th degree) before going to bed, but lowess is even better. This allows you to visualize the original distribution (its approximation) and choose the model accordingly. Thanks for the tests too. Second - I've calculated CDF incorrectly (embarassing) - thanks for pointing it out. $\endgroup$ – loard Mar 31 '14 at 7:15
2
$\begingroup$

It would be good to study a standard text on logistic regression. Among other things, you are confusing full conditioning ($Pr[Y=1 | X=x]$ as the logistic model estimates) with partial conditioning ($Pr[Y=1 | X \leq x]$). But more importantly, the model you have fitted is non-stochastic, i.e., the winner is determined with certainty from the two scores.

$\endgroup$
  • $\begingroup$ I've made an honest attempt, but I couldn't find a right example. The model does not predict with certainty (I believe) - red dots are only fitted data = prediction from the model and "score" = the position blue dots would have if forced onto regression line (forgot the technical term here). $\endgroup$ – loard Mar 30 '14 at 18:16
  • $\begingroup$ What I was trying to say is that if you know both the Ravens' and their opponents score, you know the winner with certainty. If you knew only one score you still know a lot, so the problem isn't as interesting as it might be. $\endgroup$ – Frank Harrell Mar 30 '14 at 18:35
  • $\begingroup$ You're completely right, the "problem" is very artificial. What I am trying to understand is how to visualize the distribution being fitted by logistic regression (to solve a more interesting problem, whose data I can't share). $\endgroup$ – loard Mar 30 '14 at 18:43
  • $\begingroup$ A more interesting example if you can find the data was in a paper in The American Statistician many years ago, where a binary logistic model was used to estimate the probability of a successful field goal as a function of the distance attempted, in American football. $\endgroup$ – Frank Harrell Mar 30 '14 at 19:33
  • $\begingroup$ Now I see. The logistic model is not fitting a distribution except for the Bernoulli. It is fitting the logistic function to the relationship between $x$ and $Pr[Y=1 | X=x]$. Think of it as a better way to estimate probabilities than dividing $x$ into intervals and plotting the proportion of $Y=1$ in each interval against the mean $x$ in the interval. $\endgroup$ – Frank Harrell Mar 30 '14 at 22:13
2
$\begingroup$

What the model gives is an estimate of the conditional (population) proportion of times $Y=1$ at each $x$.

If you plot the actual fitted model against the data it was fitted to:

plot(ravenWinNum~ravenScore, data=ravensData)
points(logRegRavens$fitted~ravenScore,ravensData,col=4,pch=16,cex=.8)

enter image description here

... it doesn't seem like a poor fit at all.

The two things you appear to be comparing don't seem to be comparable, but perhaps you could more clearly explain why they should be.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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