[R] ordinary polynomial coefficients from orthogonal polynomials?

Frank E Harrell Jr f.harrell at vanderbilt.edu
Wed Jun 15 05:06:35 CEST 2005


Prof Brian Ripley wrote:
> 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.

Right - I carry several digits of precision when I do this.

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

The main application I think of is when we publish fitted models, but it 
wouldn't be that bad to restate fitted orthogonal polynomials in simpler 
notation.  -Frank

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


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University




More information about the R-help mailing list