fitting a quadratic function - poly?

[Stefan Uhmann-2]

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

______________________________________________
[hidden email] 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.

[Duncan Murdoch]

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
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.

[PIKAL Petr]

Hi

[hidden email] 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
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.

[Bert Gunter]

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
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.

______________________________________________
[hidden email] 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.

[PIKAL Petr]

Hi


Bert Gunter <[hidden email]> 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 the

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

Regards
Petr

> answer is, "yes you can." See ?SafePrediction  for details. -- Bert
>
> Regards
> Petr
>
> >
> > Thanks in advance,
> > Stefan
> >
> > ______________________________________________
> > [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.

[Bert Gunter]

 

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
> >
> > ______________________________________________
> > [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.

[PIKAL Petr]

Hi

Bert Gunter <[hidden email]> 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
> > >
> > > ______________________________________________
> > > [hidden email] 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.
> >
> > ______________________________________________
> > [hidden email] 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.
> >
>
>

______________________________________________
[hidden email] 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.