[R] Confidence intervals and polynomial fits

[Ben Haller]

On May 6, 2011, at 12:31 PM, David Winsemius wrote:

> On May 6, 2011, at 11:35 AM, Ben Haller wrote:
> 
>> Hi all!  I'm getting a model fit from glm() (a binary logistic regression fit, but I don't think that's important) for a formula that contains powers of the explanatory variable up to fourth.  So the fit looks something like this (typing into mail; the actual fit code is complicated because it involves step-down and so forth):
>> 
>> x_sq <- x * x
>> x_cb <- x * x * x
>> x_qt <- x * x * x * x
>> model <- glm(outcome ~ x + x_sq + x_cb + x_qt, binomial)
> 
> It might have been easier with:
> 
> model <- glm(outcome ~ poly(x, 4) , binomial)

  Interesting.  The actual model I'm fitting has lots more terms, and needs to be able to be stepped down; sometimes some of the terms of the polynomials will get dropped while others get retained, for example.  But more to the point, poly() seems to be doing something tricky involving "orthogonal polynomials" that I don't understand.  I don't think I want whatever that is; I want my original variable x, plus its powers.  For example, if I do:

x < runif(10)
poly(x, 3)

the columns I get are not x, x^2, x^3; they are something else.  So the poly() fit is not equivalent to my fit.

> Since the authors of confint might have been expecting a poly() formulation the results might be more reliable. I'm just using lm() but I think the conclusions are more general:
> 
> ...
> 
> Coefficients:
>            Estimate Std. Error t value Pr(>|t|)
> (Intercept)  0.45737    0.52499   0.871   0.3893
> x           -0.75989    1.15080  -0.660   0.5131
> x2           1.30987    0.67330   1.945   0.0594 .
> x3          -0.03559    0.11058  -0.322   0.7494
> 
> ...
> 
> Coefficients:
>            Estimate Std. Error t value Pr(>|t|)
> (Intercept)   5.4271     0.1434  37.839  < 2e-16 ***
> poly(x, 3)1  30.0235     0.9184  32.692  < 2e-16 ***
> poly(x, 3)2   8.7823     0.9184   9.563 1.53e-11 ***
> poly(x, 3)3  -0.2956     0.9184  -0.322    0.749

  Here, in your illustration, is underscored what I mean.  Whatever these orthogonal polynomial terms are that you're using, they are clearly not the original x, x^2 and x^3, and they're giving you a different fit than those do.  I probably ought to learn about this technique, since it looks interesting; but for my purposes I need the fit to actually be in terms of x, since x is my explanatory variable.  And the fit I'm getting is highly significant (all terms < 2e-16), so the lack of fit problem you're illustrating does not seem to apply to my case.

  Or maybe I'm totally misunderstanding your point...?  :->

  Thanks!

Ben Haller
McGill University

http://biology.mcgill.ca/grad/ben/