Re: [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)36.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.80025.61414.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 -- Brian D. Ripley, [EMAIL PROTECTED] 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, UKFax: +44 1865 272595 __ 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
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
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. 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)36.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.80025.61414.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 __ 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
Re: [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)36.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.80025.61414.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 -- Roger D. Peng http://www.biostat.jhsph.edu/~rpeng/ __ 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
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
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)36.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.80025.61414.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 -- Brian D. Ripley, [EMAIL PROTECTED] 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, UKFax: +44 1865 272595 __ 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
Re: [R] ordinary polynomial coefficients from orthogonal polynomials?
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)36.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.80025.61414.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
