Confidence intervals and polynomial fits

FWIW:

Fitting higher order polynomials (say > 2) is almost always a bad idea.

See e.g.  the Hastie, Tibshirani, et. al book on "Statistical
Learning" for a detailed explanation why. The Wikipedia entry on
"smoothing splines" also contains a brief explanation, I believe.

Your ~0 P values for the coefficients also suggest problems/confusion
(!) -- perhaps you need to consider something along the lines of
"functional data analysis"  for your analysis.

Having no knowledge of your issues, these remarks are entirely
speculative and may well be wrong. So feel free to dismiss.
Nevertheless, you may find it useful to consult your local
statistician for help.

Cheers,
Bert

P.S. The raw = TRUE option of poly will fit raw, not orthogonal,
polynomials. The fitted projections will be the same up to numerical
error whichever basis is chosen, of course.

On Fri, May 6, 2011 at 10:08 AM, Ben Haller <rhelp at sticksoftware.com> wrote:

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/

-- 
"Men by nature long to get on to the ultimate truths, and will often
be impatient with elementary studies or fight shy of them. If it were
possible to reach the ultimate truths without the elementary studies
usually prefixed to them, these would not be preparatory studies but
superfluous diversions."

-- Maimonides (1135-1204)

Bert Gunter
Genentech Nonclinical Biostatistics