[R] ordinary polynomial coefficients from orthogonal polynomials?

[Prof Brian Ripley]

On Tue, 14 Jun 2005, Frank E Harrell Jr wrote:

> Prof Brian Ripley wrote:
>> On Tue, 14 Jun 2005, James Salsman wrote:
>> 
>> 
>>> How can ordinary polynomial coefficients be calculated
>>> from an orthogonal polynomial fit?
>> 
>> 
>> Why would you want to do that?  predict() is perfectly happy with an
>> orthogonal polynomial fit and the `ordinary polynomial coefficients' are 
>> rather badly determined in your example since the design matrix has a very 
>> high condition number.
>
> Brian - I don't fully see the relevance of the high condition number nowadays 
> unless the predictor has a really bad origin.  Orthogonal polynomials are a 
> mess for most people to deal with.

It means that if you write down the coeffs to a few places and then try to 
reproduce the predictions you will do badly.  The perturbation analysis 
depends on the condition number, and so is saying that the predictions are 
dependent on fine details of the coefficients.

Using (year-2000)/1000 or (year - 1970)/1000 would be a much better idea.

Why do `people' need `to deal with' these, anyway.  We have machines to do 
that.

>
> Frank
>
>> 
>> 
>>> I'm trying to do something like find a,b,c,d from
>>> lm(billions ~ a+b*decade+c*decade^2+d*decade^3)
>>> but that gives:  "Error in eval(expr, envir, enclos) :
>>> Object "a" not found"
>> 
>> 
>> You could use
>> 
>> lm(billions ~ decade + I(decade^2) + I(decade^3))
>> 
>> except that will be numerically inaccurate, since
>> 
>> 
>>> m <- model.matrix(~ decade + I(decade^2) + I(decade^3))
>>> kappa(m)
>> 
>> [1] 3.506454e+16
>> 
>> 
>> 
>> 
>>>> decade <- c(1950, 1960, 1970, 1980, 1990)
>>>> billions <- c(3.5, 5, 7.5, 13, 40)
>>>> # source: http://www.ipcc.ch/present/graphics/2001syr/large/08.17.jpg
>>>> 
>>>> pm <- lm(billions ~ poly(decade, 3))
>>>> 
>>>> plot(decade, billions, xlim=c(1950,2050), ylim=c(0,1000),
>>> 
>>> main="average yearly inflation-adjusted dollar cost of extreme weather
>>> events worldwide")
>>> 
>>>> curve(predict(pm, data.frame(decade=x)), add=TRUE)
>>>> # output: http://www.bovik.org/storms.gif
>>>> 
>>>> summary(pm)
>>> 
>>> Call:
>>> lm(formula = billions ~ poly(decade, 3))
>>> 
>>> Residuals:
>>>      1       2       3       4       5
>>> 0.2357 -0.9429  1.4143 -0.9429  0.2357
>>> 
>>> Coefficients:
>>>                 Estimate Std. Error t value Pr(>|t|)
>>> (Intercept)        13.800      0.882  15.647   0.0406 *
>>> poly(decade, 3)1   25.614      1.972  12.988   0.0489 *
>>> poly(decade, 3)2   14.432      1.972   7.318   0.0865 .
>>> poly(decade, 3)3    6.483      1.972   3.287   0.1880
>>> ---
>>> Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
>>> 
>>> Residual standard error: 1.972 on 1 degrees of freedom
>>> Multiple R-Squared: 0.9957,     Adjusted R-squared: 0.9829
>>> F-statistic: 77.68 on 3 and 1 DF,  p-value: 0.08317
>>> 
>>> 
>>>> pm
>>> 
>>> Call:
>>> lm(formula = billions ~ poly(decade, 3))
>>> 
>>> Coefficients:
>>>     (Intercept)  poly(decade, 3)1  poly(decade, 3)2  poly(decade, 3)3
>>>          13.800            25.614            14.432             6.483
>>> 
>>> ______________________________________________
>>> R-help at stat.math.ethz.ch mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide! 
>>> http://www.R-project.org/posting-guide.html
>>> 
>> 
>> 
>
>
> -- 
> Frank E Harrell Jr   Professor and Chair           School of Medicine
>                     Department of Biostatistics   Vanderbilt University
>
>

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