[R] Confidence intervals and polynomial fits

[Prof Brian Ripley]

On Fri, 6 May 2011, Bert Gunter wrote:

> 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.

That is the main piece of advice I would have given.  But if you must 
DIY, consider the merits of orthogonal polynomials.  Computing 
individual confidence intervals for highly correlated coefficients is 
very dubious practice.  Without the example the posting guide asked 
for, we can only guess if that is what is happening.

>
> 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/
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
>
> -- 
> "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
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595