5
$\begingroup$

This is a follow-up question from this post, here: https://stackoverflow.com/questions/14423325/confidence-intervals-for-predictions-from-logistic-regression

The answer from @Gavin is excellent, but I have some additional questions which I think would be useful for others. I am working with a Poisson model, so basically it is the same approach described in the other post, only family=poisson instead of family=binomial.

To my first question: @Gavin writes:

mod <- glm(y ~ x, data = foo, family = binomial)
preddat <- with(foo, data.frame(x = seq(min(x), max(x), length = 100))
preds <- predict(mod, newdata = preddata, type = "link", se.fit = TRUE)

What is the point of the second line there? Is it necessary to create a data.frame with minimum and maximum of the explanatory variable? Could I not, for some explanatory variable(s) x (stored in some data frame data), just go from the first line and directly to the third?

To my second question: In the beginning of his answer @Gavin writes:

The usual way is to compute a confidence interval on the scale of the linear predictor, where things will be more normal (Gaussian) and then apply the inverse of the link function to map the confidence interval from the linear predictor scale to the response scale.

Why are "things" more normal on the scale of the linear predictor(s)? Is this also the case when I do my Poisson regression? I assume the reason for using critical value 1.96 when constructing the CI's, is because of the assumptions that "things" are normal. Can somebody explain this further?

My third question:

Is there a relationship between the standard deviation which we get by using se.fit=TRUE in predict() and the standard deviations of the coefficients of the explanatory variables, which we simply get from summary(mod)? (mod is some glm object)

$\endgroup$

migrated from stackoverflow.com May 1 '14 at 16:37

This question came from our site for professional and enthusiast programmers.

  • $\begingroup$ This question strikes me as on the border between here (SO) & stats.SE (Cross Validated). In particular, the 2nd Q in this post is definitely statistical. $\endgroup$ – gung Apr 24 '14 at 15:27
  • $\begingroup$ I know, I was very uncertain whether to put it here or stats.SE. I im truly sorry if I was wrong. $\endgroup$ – Erosennin Apr 24 '14 at 15:33
  • $\begingroup$ You aren't necessarily wrong, @ErosRam, & there's no need to apologize. I just wanted to make the comment. $\endgroup$ – gung Apr 24 '14 at 15:38
5
$\begingroup$

To address Q1, lets start by making some data to play with:

lo.to.p <- function(lo){  # this function will convert log odds to probabilities
  o <- exp(lo)            # we get odds by exponentiating log odds
  p <- o/(o+1)            # we convert to probabilities
  return(p)
}

set.seed(90)                        # this makes the example reproducible
x   <- runif(100, min=0, max=100)   # I generate some x data from a uniform dist
lo  <- -.5 + .1*x                   # this is the linear predictor
p   <- lo.to.p(lo)                  # converting log odds to probabilities
y   <- rbinom(100, size=1, prob=p)  # generating observed y values
foo <- data.frame(x=x, y=y)

  # @Gavin's code:
mod     <- glm(y ~ x, data=foo, family=binomial)
preddat <- with(foo, data.frame(x=seq(min(x), max(x), length=100)))
preds   <- predict(mod, newdata=preddat, type="link", se.fit=TRUE)

Now, why not try to get predicted values and a confidence interval / band by just using the original data:

preds2  <- predict(mod, newdata=foo$x, type="link", se.fit=TRUE)

That throws an error, because predict() needs the newdata argument to get a data frame:

# Error in eval(predvars, data, env) : 
#   numeric 'envir' arg not of length one

So let's try with the original data as a data frame:

preds3  <- predict(mod, newdata=data.frame(x=foo$x), type="link", se.fit=TRUE)

That time it worked, so let's see what the output looks like (I used our lo.to.p() function to convert the output from predict to predicted probabilities as @Gavin suggested, note that you can also use predict with type="response" to do that automatically):

enter image description here

Using the original data frame yields a garbled mess. You can sort the data first, which works OK in this case, but generally is not as smooth / pretty. To better show the effect of this strategy, I slightly augmented the data and model. Here's the code for the sorted version:

foo2    <- with(foo, data.frame(x=c(x, -100), y=c(y,0)))
mod2    <- glm(y~x, data=foo2, family=binomial)
preds4  <- predict(mod2, newdata=data.frame(x=sort(foo2$x)), type="link", 
                   se.fit=TRUE)

Regarding Q2, the statistical theory behind generalized linear models (GLiMs) assumes that the sampling distribution of a parameter estimate is asymptotically normally distributed (i.e., 'at infinity'). It is well known that this is not necessarily true for small samples, but the sampling distribution may be 'normal enough'. At any rate, this is (possibly) true on the scale of the linear predictor, which I call lo above; but the link function is a non-linear transformation, it isn't necessarily true on the response scale. To use an easy example, the normal distribution goes to infinity on both sides, but the response scale is bounded at 0 and 1. Moreover, all of these points hold for the Poisson distribution just like the binomial. Although it's not exactly the same topic, it may help to read my answer here: difference between logit and probit models because it provides a lot of information about link functions and GLiMs that may help with the larger conceptual framework.

For Q3, yes there is a relationship between the SEs of your coefficients and the width confidence band, but the confidence band is a little more complicated. The width of the confidence band grows as you move left or right away from the mean of x. (You can get the general idea from my answer here: linear regression prediction interval.) On the other hand, with a GLiM, the width of the confidence band also depends on the predicted value. To more easily see these effects, we can look at the confidence band for our original model on the scale of the linear predictor, and for a second model where there is no effect of x. Here's the second model:

y2      <- rbinom(100, size=1, prob=.5)
mod2    <- glm(y2~x, family=binomial)
preds5  <- predict(mod2, newdata=data.frame(x=sort(foo$x)), type="link", 
                   se.fit=TRUE)

Here's what they look like:

enter image description here

$\endgroup$
  • $\begingroup$ Awesome answer, thank you :) With regard to Q1: so we sort the values (or arrange them in a sequence) to get more smooth predictions - but if I am just interested in the predictions themselves, in no particular order, I do not have to arrange them, correct? $\endgroup$ – Erosennin Apr 25 '14 at 8:21
  • $\begingroup$ @ErosRam, you need Gavin's 2nd line of code to make a plot. I tweaked the code / plot to better demonstrate this. If you don't want to make a plot, but get the predicted probabilities (or lambda for Poisson reg) +/- SEs for your data, then you would use your original data & the sorting is irrelevant. $\endgroup$ – gung Apr 25 '14 at 17:07
  • 1
    $\begingroup$ To get predictions for factors, you use the same formula (at least for linear models), or, more likely a multidimensional version of it in matrix form. I don't remember the exact formula off the top of my head, but these are standard in textbooks. Moreover you would need a Poisson or logistic (etc) specific version, b/c the variance scales w/ the predicted value (note my answer to Q3 above). You might just ask for the specific formula you need on CV. $\endgroup$ – gung Apr 25 '14 at 17:08

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.