[R] Orthogonal polynomials used by R
Ashim Kapoor
@@h|mk@poor @end|ng |rom gm@||@com
Thu Nov 28 07:16:42 CET 2019
Dear Bert,
OK and thank you.
@Fox, John <jfox using mcmaster.ca> will be grateful for an offline reply.
Best,
Ashim
On Thu, Nov 28, 2019 at 11:43 AM Bert Gunter <bgunter.4567 using gmail.com> wrote:
> Statistical questions are generally off topic on this list. Try
> stats.stackexchange.com instead.
>
> But FWIW, I recommend that you work with someone with expertise in time
> series analysis, as your efforts to shake and bake your own methodology
> seem rather unwise to me.
>
> Cheers,
> Bert
>
> Bert Gunter
>
> "The trouble with having an open mind is that people keep coming along and
> sticking things into it."
> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
>
>
> On Wed, Nov 27, 2019 at 9:47 PM Ashim Kapoor <ashimkapoor using gmail.com>
> wrote:
>
>> Dear Peter and John,
>>
>> Many thanks for your prompt replies.
>>
>> Here is what I was trying to do. I was trying to build a statistical
>> model
>> of a given time series using Box Jenkins methodology. The series has 93
>> data points. Before I analyse the ACF and PACF, I am required to de-trend
>> the series. The series seems to have an upward trend. I wanted to find out
>> what order polynomial should I fit the series
>> without overfitting. For this I want to use orthogonal polynomials(I
>> think
>> someone on the internet was talking about preventing overfitting by using
>> orthogonal polynomials) . This seems to me as a poor man's cross
>> validation.
>>
>> So my plan is to keep increasing the degree of the orthogonal polynomials
>> till the coefficient of the last orthogonal polynomial becomes
>> insignificant.
>>
>> Note : If I do NOT use orthogonal polynomials, I will overfit the data set
>> and I don't think that is a good way to detect the true order of the
>> polynomial.
>>
>> Also now that I have detrended the series and built an ARIMA model of the
>> residuals, now I want to forecast. For this I need to use the original
>> polynomials and their coefficients.
>>
>> I hope I was clear and that my methodology is ok.
>>
>> I have another query here :-
>>
>> Note : If I used cross-validation to determine the order of the
>> polynomial,
>> I don't get a clear answer.
>>
>> See here :-
>> library(boot)
>> mydf = data.frame(cbind(gdp,x))
>> d<-(c(
>> cv.glm(data = mydf,glm(gdp~x),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~poly(x,2)),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~poly(x,3)),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~poly(x,4)),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~poly(x,5)),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~poly(x,6)),K=10)$delta[1]))
>> print(d)
>> ## [1] 2.178574e+13 7.303031e+11 5.994783e+11 4.943586e+11 4.596648e+11
>> ## [6] 4.980159e+11
>>
>> # Here it chooses 5. (but 4 and 5 are kind of similar).
>>
>>
>> d1 <- (c(
>> cv.glm(data = mydf,glm(gdp~1+x),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~1+x+x^2),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~1+x+x^2+x^3),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~1+x+x^2+x^3+x^4),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~1+x+x^2+x^3+x^4+x^5),K=10)$delta[1],
>> cv.glm(data = mydf,glm(gdp~1+x+x^2+x^3+x^4+x^5+x^6),K=10)$delta[1]))
>>
>> print(d1)
>> ## [1] 2.149647e+13 2.253999e+13 2.182175e+13 2.177170e+13 2.198675e+13
>> ## [6] 2.145754e+13
>>
>> # here it chooses 1 or 6
>>
>> Query : Why does it choose 1? Notice : Is this just round off noise /
>> noise
>> due to sampling error created by Cross Validation when it creates the K
>> folds? Is this due to the ill conditioned model matrix?
>>
>> Best Regards,
>> Ashim.
>>
>>
>>
>>
>>
>> On Wed, Nov 27, 2019 at 10:41 PM Fox, John <jfox using mcmaster.ca> wrote:
>>
>> > Dear Ashim,
>> >
>> > Orthogonal polynomials are used because they tend to produce more
>> accurate
>> > numerical computations, not because their coefficients are
>> interpretable,
>> > so I wonder why you're interested in the coefficients.
>> >
>> > The regressors produced are orthogonal to the constant regressor and are
>> > orthogonal to each other (and in fact are orthonormal), as it's simple
>> to
>> > demonstrate:
>> >
>> > ------- snip -------
>> >
>> > > x <- 1:93
>> > > y <- 1 + x + x^2 + x^3 + x^4 + rnorm(93)
>> > > (m <- lm(y ~ poly(x, 4)))
>> >
>> > Call:
>> > lm(formula = y ~ poly(x, 4))
>> >
>> > Coefficients:
>> > (Intercept) poly(x, 4)1 poly(x, 4)2 poly(x, 4)3 poly(x, 4)4
>> > 15574516 172715069 94769949 27683528 3429259
>> >
>> > > X <- model.matrix(m)
>> > > head(X)
>> > (Intercept) poly(x, 4)1 poly(x, 4)2 poly(x, 4)3 poly(x, 4)4
>> > 1 1 -0.1776843 0.2245083 -0.2572066 0.27935949
>> > 2 1 -0.1738216 0.2098665 -0.2236579 0.21862917
>> > 3 1 -0.1699589 0.1955464 -0.1919525 0.16390514
>> > 4 1 -0.1660962 0.1815482 -0.1620496 0.11487597
>> > 5 1 -0.1622335 0.1678717 -0.1339080 0.07123722
>> > 6 1 -0.1583708 0.1545171 -0.1074869 0.03269145
>> >
>> > > zapsmall(crossprod(X))# X'X
>> > (Intercept) poly(x, 4)1 poly(x, 4)2 poly(x, 4)3 poly(x, 4)4
>> > (Intercept) 93 0 0 0 0
>> > poly(x, 4)1 0 1 0 0 0
>> > poly(x, 4)2 0 0 1 0 0
>> > poly(x, 4)3 0 0 0 1 0
>> > poly(x, 4)4 0 0 0 0 1
>> >
>> > ------- snip -------
>> >
>> > If for some not immediately obvious reason you're interested in the
>> > regression coefficients, why not just use a "raw" polynomial:
>> >
>> > ------- snip -------
>> >
>> > > (m1 <- lm(y ~ poly(x, 4, raw=TRUE)))
>> >
>> > Call:
>> > lm(formula = y ~ poly(x, 4, raw = TRUE))
>> >
>> > Coefficients:
>> > (Intercept) poly(x, 4, raw = TRUE)1 poly(x, 4, raw =
>> TRUE)2
>> > poly(x, 4, raw = TRUE)3
>> > 1.5640 0.8985
>> 1.0037
>> > 1.0000
>> > poly(x, 4, raw = TRUE)4
>> > 1.0000
>> >
>> > ------- snip -------
>> >
>> > These coefficients are simply interpretable but the model matrix is more
>> > poorly conditioned:
>> >
>> > ------- snip -------
>> >
>> > > head(X1)
>> > (Intercept) poly(x, 4, raw = TRUE)1 poly(x, 4, raw = TRUE)2 poly(x, 4,
>> > raw = TRUE)3
>> > 1 1 1 1
>> > 1
>> > 2 1 2 4
>> > 8
>> > 3 1 3 9
>> > 27
>> > 4 1 4 16
>> > 64
>> > 5 1 5 25
>> > 125
>> > 6 1 6 36
>> > 216
>> > poly(x, 4, raw = TRUE)4
>> > 1 1
>> > 2 16
>> > 3 81
>> > 4 256
>> > 5 625
>> > 6 1296
>> > > round(cor(X1[, -1]), 2)
>> > poly(x, 4, raw = TRUE)1 poly(x, 4, raw = TRUE)2
>> > poly(x, 4, raw = TRUE)3
>> > poly(x, 4, raw = TRUE)1 1.00 0.97
>> > 0.92
>> > poly(x, 4, raw = TRUE)2 0.97 1.00
>> > 0.99
>> > poly(x, 4, raw = TRUE)3 0.92 0.99
>> > 1.00
>> > poly(x, 4, raw = TRUE)4 0.87 0.96
>> > 0.99
>> > poly(x, 4, raw = TRUE)4
>> > poly(x, 4, raw = TRUE)1 0.87
>> > poly(x, 4, raw = TRUE)2 0.96
>> > poly(x, 4, raw = TRUE)3 0.99
>> > poly(x, 4, raw = TRUE)4 1.00
>> >
>> > ------- snip -------
>> >
>> > The two parametrizations are equivalent, however, in that they represent
>> > the same regression surface, and so, e.g., produce the same fitted
>> values:
>> >
>> > ------- snip -------
>> >
>> > > all.equal(fitted(m), fitted(m1))
>> > [1] TRUE
>> >
>> > ------- snip -------
>> >
>> > Because one is usually not interested in the individual coefficients of
>> a
>> > polynomial there usually isn't a reason to prefer one parametrization to
>> > the other on the grounds of interpretability, so why do you need to
>> > interpret the regression equation?
>> >
>> > I hope this helps,
>> > John
>> >
>> > -----------------------------
>> > John Fox, Professor Emeritus
>> > McMaster University
>> > Hamilton, Ontario, Canada
>> > Web: http::/socserv.mcmaster.ca/jfox
>> >
>> > > On Nov 27, 2019, at 10:17 AM, Ashim Kapoor <ashimkapoor using gmail.com>
>> > wrote:
>> > >
>> > > Dear Petr,
>> > >
>> > > Many thanks for the quick response.
>> > >
>> > > I also read this:-
>> > > https://en.wikipedia.org/wiki/Discrete_orthogonal_polynomials
>> > >
>> > > Also I read in ?poly:-
>> > > The orthogonal polynomial is summarized by the coefficients, which
>> > > can be used to evaluate it via the three-term recursion given in
>> > > Kennedy & Gentle (1980, pp. 343-4), and used in the ‘predict’ part
>> > > of the code.
>> > >
>> > > I don't have access to the mentioned book.
>> > >
>> > > Out of curiosity, what is the name of the discrete orthogonal
>> polynomial
>> > > used by R ?
>> > > What discrete measure is it orthogonal with respect to ?
>> > >
>> > > Many thanks,
>> > > Ashim
>> > >
>> > >
>> > >
>> > >
>> > > On Wed, Nov 27, 2019 at 6:11 PM PIKAL Petr <petr.pikal using precheza.cz>
>> > wrote:
>> > >
>> > >> You could get answer quickly by searching net.
>> > >>
>> > >>
>> > >>
>> >
>> https://stackoverflow.com/questions/39031172/how-poly-generates-orthogonal-p
>> > >> olynomials-how-to-understand-the-coefs-ret/39051154#39051154
>> > >> <
>> >
>> https://stackoverflow.com/questions/39031172/how-poly-generates-orthogonal-polynomials-how-to-understand-the-coefs-ret/39051154#39051154
>> > >
>> > >>
>> > >> Cheers
>> > >> Petr
>> > >>
>> > >>> -----Original Message-----
>> > >>> From: R-help <r-help-bounces using r-project.org> On Behalf Of Ashim
>> Kapoor
>> > >>> Sent: Wednesday, November 27, 2019 12:55 PM
>> > >>> To: R Help <r-help using r-project.org>
>> > >>> Subject: [R] Orthogonal polynomials used by R
>> > >>>
>> > >>> Dear All,
>> > >>>
>> > >>> I have created a time trend by doing x<-1:93 because I have a time
>> > series
>> > >>> with 93 data points. Next I did :-
>> > >>>
>> > >>> y = lm(series ~ poly(x,4))$residuals
>> > >>>
>> > >>> to detrend series.
>> > >>>
>> > >>> I choose this 4 as the order of my polynomial using cross
>> validation/
>> > >>> checking the absence of trend in the residuals so I think I have not
>> > >> overfit
>> > >>> this series.
>> > >>>
>> > >>> I wish to document the formula of poly(x,4). I am not able to find
>> it
>> > in
>> > >> ?poly
>> > >>>
>> > >>> Can someone please tell me what the formula for the orthogonal
>> > >>> polynomial used by R is ?
>> > >>>
>> > >>> Thank you,
>> > >>> Ashim
>> > >>>
>> > >>> [[alternative HTML version deleted]]
>> > >>>
>> > >>> ______________________________________________
>> > >>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> > >>> 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.
>> > >>
>> > >
>> > > [[alternative HTML version deleted]]
>> > >
>> > > ______________________________________________
>> > > R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> > > 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.
>> >
>> >
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>>
>
[[alternative HTML version deleted]]
More information about the R-help
mailing list