Hi List, I can not get my head around the following problem. I want to fit a quadratic function to some data and stumbled across poly(). What exactly does it, i.e. why are there different results for fit1 and fit2? x = seq(-10, 10) y = x^2 fit1 = lm(y ~ x + I(x^2)) fit2 = lm(y ~ poly(x, 2)) plot(x,y) lines(x, fit1$fitted.values, col = 2) lines(x, fit2$fitted.values, col = 3) round(fit1$coefficients, 2) round(fit2$coefficients, 2) Thanks in advance, Stefan
On 14/04/2010 11:12 AM, Stefan Uhmann wrote:> Hi List, > > I can not get my head around the following problem. I want to fit a > quadratic function to some data and stumbled across poly(). What exactly > does it, i.e. why are there different results for fit1 and fit2? > > x = seq(-10, 10) > y = x^2 > > fit1 = lm(y ~ x + I(x^2)) > fit2 = lm(y ~ poly(x, 2)) > > plot(x,y) > lines(x, fit1$fitted.values, col = 2) > lines(x, fit2$fitted.values, col = 3) >These look the same to me.> round(fit1$coefficients, 2) > round(fit2$coefficients, 2) >These look different, because poly uses orthogonal polynomials, a different parametrization. You can see the difference if you ask for model.matrix(fit1) and model.matrix(fit2). (You can plot these using matplot(model.matrix(fit1)), etc.) Duncan Murdoch> Thanks in advance, > Stefan > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >
Hi r-help-bounces at r-project.org napsal dne 14.04.2010 17:12:51:> Hi List, > > I can not get my head around the following problem. I want to fit a > quadratic function to some data and stumbled across poly(). What exactly> does it, i.e. why are there different results for fit1 and fit2? > > x = seq(-10, 10) > y = x^2 > > fit1 = lm(y ~ x + I(x^2)) > fit2 = lm(y ~ poly(x, 2)) > > plot(x,y) > lines(x, fit1$fitted.values, col = 2) > lines(x, fit2$fitted.values, col = 3)> round(fitted(fit1)-fitted(fit2),5)1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 results are same.> > round(fit1$coefficients, 2) > round(fit2$coefficients, 2)Coefficients are different as you fit different values. See ?poly poly(-10:10,2) I believe that others give you better explanation. So you can not use coefficients evaluated by lm(.~poly(...)) directly. Regards Petr> > Thanks in advance, > Stefan > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code.
Below. -- Bert Bert Gunter Genentech Nonclinical Statistics Coefficients are different as you fit different values. See ?poly poly(-10:10,2) I believe that others give you better explanation. So you can not use coefficients evaluated by lm(.~poly(...)) directly. -- Well, it depends what you mean by "use...directly." But I think the answer is, "yes you can." See ?SafePrediction for details. -- Bert Regards Petr> > Thanks in advance, > Stefan > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code.______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Hi Bert Gunter <gunter.berton at gene.com> napsal dne 14.04.2010 18:01:52:> Below. > > -- Bert > > > Bert Gunter > Genentech Nonclinical Statistics > > > Coefficients are different as you fit different values. See > > ?poly > > poly(-10:10,2) > > I believe that others give you better explanation. So you can not use > coefficients evaluated by lm(.~poly(...)) directly. > > -- Well, it depends what you mean by "use...directly." But I think theI mean that you can use fit<- lm(y~x+I(x^2)) coef(fit)[1] + coef(fit)[2]*x + coef(fit)[3]*x^2 but you can not use fit<- lm(y~poly(x,2)) coef(fit)[1] + coef(fit)[2]*x + coef(fit)[3]*x^2 to compute y. Regards Petr> answer is, "yes you can." See ?SafePrediction for details. -- Bert > > Regards > Petr > > > > > Thanks in advance, > > Stefan > > > > ______________________________________________ > > R-help at r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. >
Petr Pikal wrote: "... I mean that you can use fit<- lm(y~x+I(x^2)) coef(fit)[1] + coef(fit)[2]*x + coef(fit)[3]*x^2 but you can not use fit<- lm(y~poly(x,2)) coef(fit)[1] + coef(fit)[2]*x + coef(fit)[3]*x^2 " (to get the fits for any x vector) -- But you **can** use ypred <- predict(fit,data.frame(x = x)) -- in **both** cases. Which is, I think, how it should be done. Cheers, Bert Bert Gunter Genentech Nonclinical Statistics> answer is, "yes you can." See ?SafePrediction for details. -- Bert > > Regards > Petr > > > > > Thanks in advance, > > Stefan > > > > ______________________________________________ > > R-help at r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. >
Hi Bert Gunter <gunter.berton at gene.com> napsal dne 14.04.2010 18:54:37:> > > Petr Pikal wrote: > > "... > > I mean that you can use > > fit<- lm(y~x+I(x^2)) > coef(fit)[1] + coef(fit)[2]*x + coef(fit)[3]*x^2 > > but you can not use > > fit<- lm(y~poly(x,2)) > coef(fit)[1] + coef(fit)[2]*x + coef(fit)[3]*x^2 " > > (to get the fits for any x vector) > > -- But you **can** use > > ypred <- predict(fit,data.frame(x = x)) > > -- in **both** cases. Which is, I think, how it should be done.I completely agree with you. However sometimes it is necessary to put an equation and estimated coefficients to a report - for that you can not use coefficients from poly() estimate directly. (Please do not beat me, I know that polynomial model is ***rarely*** physically viable and personally I never use it until it it has not a physical sense - like free fall) Regards Petr> > Cheers, > Bert > > Bert Gunter > Genentech Nonclinical Statistics > > > > > answer is, "yes you can." See ?SafePrediction for details. -- Bert > > > > Regards > > Petr > > > > > > > > Thanks in advance, > > > Stefan > > > > > > ______________________________________________ > > > R-help at r-project.org mailing list > > > https://stat.ethz.ch/mailman/listinfo/r-help > > > PLEASE do read the posting guide > > http://www.R-project.org/posting-guide.html > > > and provide commented, minimal, self-contained, reproducible code. > > > > ______________________________________________ > > R-help at r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > > >