James Salsman
2005-Jun-14 09:03 UTC
[R] ordinary polynomial coefficients from orthogonal polynomials?
How can ordinary polynomial coefficients be calculated
from an orthogonal polynomial fit?
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"
> 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
Prof Brian Ripley
2005-Jun-14 10:35 UTC
[R] ordinary polynomial coefficients from orthogonal polynomials?
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.> 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 >-- 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
Roger D. Peng
2005-Jun-14 13:48 UTC
[R] ordinary polynomial coefficients from orthogonal polynomials?
I think this is covered in the FAQ: http://cran.r-project.org/doc/FAQ/R-FAQ.html#Models -roger James Salsman wrote:> How can ordinary polynomial coefficients be calculated > from an orthogonal polynomial fit? > > 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" > > > 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 >-- Roger D. Peng http://www.biostat.jhsph.edu/~rpeng/