From smckinney at bccrc.ca  Fri Apr 25 05:26:55 2008
From: smckinney at bccrc.ca (Steven McKinney)
Date: Thu, 24 Apr 2008 20:26:55 -0700
Subject: [Rd] poly() can exceed degree k - 1 for k distinct points
	(PR#11251)
References: <mailman.22.1209031205.8468.r-devel@r-project.org>
	<FEFDBEC5-CF9D-465E-A075-EB071252B66D@anu.edu.au>
Message-ID: <0BE438149FF2254DB4199E2682C8DFEB0328A281@crcmail1.BCCRC.CA>


I haven't parsed the source to fully understand
the 'normalization constants' returned in the
poly output component 
  attr(,"coefs")$norm2
but notice that the first 6 are non-zero and
the last 8 are smaller than machine precision.
Some kind of useful warning to the user could
ensue upon evaluation of such normalization
constants?



> x = rep(1:5,3)
> x
 [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
> poly(x, 12)
                  1          2             3           4           5           6           7           8           9          10
 [1,] -3.651484e-01  0.3086067 -1.825742e-01  0.06900656  0.50410540 -0.04690012  0.16698778  0.03700268 -0.21853454 -0.50410540
 [2,] -1.825742e-01 -0.1543033  3.651484e-01 -0.27602622  0.11156408 -0.48312338 -0.10874588 -0.02250863  0.27412822 -0.11156408
 [3,] -6.437201e-18 -0.3086067 -1.301129e-16  0.41403934 -0.17911692 -0.14493516  0.51592172  0.11434917  0.03151684  0.17911692
 [4,]  1.825742e-01 -0.1543033 -3.651484e-01 -0.27602622  0.07326369  0.11056157 -0.07141299  0.52748012  0.18001893 -0.07326369
 [5,]  3.651484e-01  0.3086067  1.825742e-01  0.06900656 -0.17121360 -0.25395954 -0.14949152  0.20036589 -0.42074883  0.17121360
 [6,] -3.651484e-01  0.3086067 -1.825742e-01  0.06900656 -0.75205270  0.02345006 -0.08349389 -0.01850134  0.10926727 -0.24794730
 [7,] -1.825742e-01 -0.1543033  3.651484e-01 -0.27602622 -0.05578204  0.74156169  0.05437294  0.01125431 -0.13706411  0.05578204
 [8,] -6.437201e-18 -0.3086067 -8.788409e-17  0.41403934  0.08955846  0.07246758 -0.75796086 -0.05717459 -0.01575842 -0.08955846
 [9,]  1.825742e-01 -0.1543033 -3.651484e-01 -0.27602622 -0.03663185 -0.05528079  0.03570650 -0.76374006 -0.09000946  0.03663185
[10,]  3.651484e-01  0.3086067  1.825742e-01  0.06900656  0.08560680  0.12697977  0.07474576 -0.10018295  0.71037442 -0.08560680
[11,] -3.651484e-01  0.3086067 -1.825742e-01  0.06900656  0.24794730  0.02345006 -0.08349389 -0.01850134  0.10926727  0.75205270
[12,] -1.825742e-01 -0.1543033  3.651484e-01 -0.27602622 -0.05578204 -0.25843831  0.05437294  0.01125431 -0.13706411  0.05578204
[13,] -6.437201e-18 -0.3086067 -8.788409e-17  0.41403934  0.08955846  0.07246758  0.24203914 -0.05717459 -0.01575842 -0.08955846
[14,]  1.825742e-01 -0.1543033 -3.651484e-01 -0.27602622 -0.03663185 -0.05528079  0.03570650  0.23625994 -0.09000946  0.03663185
[15,]  3.651484e-01  0.3086067  1.825742e-01  0.06900656  0.08560680  0.12697977  0.07474576 -0.10018295 -0.28962558 -0.08560680
               11          12
 [1,]  0.04690012  0.16698778
 [2,]  0.48312338 -0.10874588
 [3,]  0.14493516  0.51592172
 [4,] -0.11056157 -0.07141299
 [5,]  0.25395954 -0.14949152
 [6,] -0.02345006 -0.08349389
 [7,]  0.25843831  0.05437294
 [8,] -0.07246758  0.24203914
 [9,]  0.05528079  0.03570650
[10,] -0.12697977  0.07474576
[11,] -0.02345006 -0.08349389
[12,] -0.74156169  0.05437294
[13,] -0.07246758 -0.75796086
[14,]  0.05528079  0.03570650
[15,] -0.12697977  0.07474576
attr(,"degree")
 [1]  1  2  3  4  5  6  7  8  9 10 11 12
attr(,"coefs")
attr(,"coefs")$alpha
 [1] 3.000000 3.000000 3.000000 3.000000 3.000000 1.314957 2.355111 2.973300 4.032925 4.323708 1.314957 2.355111

attr(,"coefs")$norm2
 [1] 1.000000e+00 1.500000e+01 3.000000e+01 4.200000e+01 4.320000e+01 2.468571e+01 1.128824e-28 1.344243e-28 7.418009e-29 9.824087e-27
[11] 1.017082e-32 2.357175e-51 9.543887e-57 6.351670e-57

attr(,"class")
[1] "poly"   "matrix"


Steven McKinney



-----Original Message-----
From: r-devel-bounces at r-project.org on behalf of John Maindonald
Sent: Thu 4/24/2008 5:24 PM
To: r-devel at r-project.org
Subject: Re: [Rd] R-devel Digest, Vol 62, Issue 24
 
Actually, this may be a useful feature!  It allows calculation of a  
basis for the orthogonal complement of the space spanned by  
model.matrix(lm(y ~ poly(x,12)).  However, the default ought surely to  
be to disallow df > k-1 in poly(x,df), where k = length(unique(x)).

John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.


On 24 Apr 2008, at 8:00 PM, r-devel-request at r-project.org wrote:

> From: russell-lenth at uiowa.edu
> Date: 24 April 2008 3:05:28 AM
> To: r-devel at stat.math.ethz.ch
> Cc: R-bugs at biostat.ku.dk
> Subject: [Rd] poly() can exceed degree k - 1 for k distinct points  
> (PR#11251)
>
>
> The poly() function can create more variables than can be fitted when
> there are replicated values.  In the example below, 'x' has only 5
> distinct values, but I can apparently fit a 12th-degree polynomial  
> with
> no error messages or even nonzero coefficients:
>
> R> x = rep(1:5,3)
> R> y = rnorm(15)
> R> lm(y ~ poly(x, 12))
>
> Call:
> lm(formula = y ~ poly(x, 12))
>
> Coefficients:
>   (Intercept)   poly(x, 12)1   poly(x, 12)2   poly(x, 12)3
>      -0.27442        0.35822       -0.26412        2.11780
>  poly(x, 12)4   poly(x, 12)5   poly(x, 12)6   poly(x, 12)7
>       1.83117       -0.09260       -0.48572        1.94030
>  poly(x, 12)8   poly(x, 12)9  poly(x, 12)10  poly(x, 12)11
>      -0.88297       -1.04556        0.74289       -0.01422
> poly(x, 12)12
>      -0.46548
>
<snip>
<snip>

> [I thought I submitted this via the website yesterday, but I can  
> find no
> trace of it.  I apologize if this is a duplicate, but I don't think  
> it is.]
> -- 
> Russell V. Lenth, Professor
> Department of Statistics
>   & Actuarial Science            (319)335-0814    FAX (319)335-3017
> The University of Iowa           russell-lenth at uiowa.edu
> Iowa City, IA 52242  USA         http://www.stat.uiowa.edu/~rlenth/


	[[alternative HTML version deleted]]

______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel



[Rd] poly() can exceed degree k - 1 for k distinct points (PR#11251)
Steven McKinney