From: Frank E Harrell Jr <f.harrell_at_vanderbilt.edu>

Date: Tue 14 Jun 2005 - 23:46:50 EST

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

>>m <- model.matrix(~ decade + I(decade^2) + I(decade^3))

*>>kappa(m)
*

*>>>
*

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

*>>
*

>

>

Date: Tue 14 Jun 2005 - 23:46:50 EST

> On Tue, 14 Jun 2005, James Salsman wrote: > >

>>How can ordinary polynomial coefficients be calculated

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

Frank

> >

>>I'm trying to do something like find a,b,c,d from

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

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

>

>

-- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University ______________________________________________ R-help@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.htmlReceived on Tue Jun 14 23:48:45 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:32:34 EST
*