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