Re: [R] RDA and trend surface regression
Dear All, PCNM and more general MEM can be found in the package spacemakeR available here (http://biomserv.univ-lyon1.fr/~dray/). As said by Jari, I do not want to put the package on CRAN. I have no no problems to include it in another package (vegan or spdep). It has a vignette which can be interesting for spdep. However, I have some problems with this vector approach for selection of variables in multivariate analysis. spdep has some function to select these eigenvectors in univariate regression, we have to find some rules for selecting in the multivariate case. I use AIC for the paper but it is quite sure that it is not correct (see comments of Jari in the vegan help page). multispati is a spatially constrained multivariate analysis. I think that this approach is more clean in a 'statistical sense' but quite different: we do not have variance but autocorrelation and so negative eigenvalues in multivariate analysis ! I have see that Roger has send an announce for a spatial-R meeting in R-Geo List. I am interested to participate to the Bergen session (I love rain ;-) ). It could be nice to speak about this subject in Bergen. Cheers, Jari Oksanen wrote: On 27 Feb 2007, at 20:55, Gavin Simpson wrote: On Tue, 2007-02-27 at 13:13 -0500, Kuhn, Max wrote: Helene, My point was only that RDA may fit a quadratic model for the terms specified in your model. The terms that you had specified were already higher order polynomials (some cubic). So a QDA classifier with the model terms that you specified my be a fifth order polynomial in the original data. I don't know the reference you cite or even the subject-matter specifics. I'm just a simple cave man (for you SNL fans). But I do know that there are more reliable ways to get nonlinear classification boundaries than using x^5. I doubt that Helene is trying to do a classification - unless you consider classification to mean that all rows/samples are in different groups (i.e. n samples therefore n groups) - which is how RDA (Redundancy Analysis) is used in ecology. You could take a look at multispati in package ade4 for a different way to handle spatial constraints. There is also the principle coordinates analysis of neighbour matrices (PCNM) method - not sure this is coded anywhere in R yet though. Here are two references that may be useful: Stéphane Dray has R code for finding PCNM matrices. Google for his name: it's not that common. I also have a copy of his function and can send it if really needed, but it may be better to check Dray's page first. Stéphane Dray says think that not all functions need be in CRAN. May be true, but I think it might help many people. There are at least three reasons why not use polynomial constraints in RDA. Max Kuhn mentioned one: polynomials typically flip wildly at margins (or they are unstable in more neutral speech). Second reason is that they are almost impossible to interpret in ordination display. The third reason is that RDA (or CCA) avoid some ordination artefacts (curvature, horseshoe, arc etc.) just because the constraints are linear: allowing them to be curved allows curved solutions. These arguments are not necessarily valid if you only want to have variance partitioning, or if you use polynomial conditions (partial out polynomial effects in Canoco language). In that case it may make sense to use quadratic (or polynomial) constraints or conditions. cheers, Jari Oksanen __ 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 and provide commented, minimal, self-contained, reproducible code. -- Stéphane DRAY ([EMAIL PROTECTED] ) Laboratoire BBE-CNRS-UMR-5558, Univ. C. Bernard - Lyon I 43, Bd du 11 Novembre 1918, 69622 Villeurbanne Cedex, France Tel: 33 4 72 43 27 57 Fax: 33 4 72 43 13 88 http://biomserv.univ-lyon1.fr/~dray/ __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] RDA and trend surface regression
Thanks a lot for your answers, I am concerned by your advice not to use polynomial constraints, or to use QDA instead of RDA. My final goal is to perform variation partitioning using partial RDA to assess the relative importance of environmental vs spatial variables. For the spatial analyses, trend surface analysis (polynomial constraints) is recommended in Legendre and Legendre 1998 (p739). Is there a better method to integrate space as an explanatory variable in a variation partitioning analyses? Also, I don't understand this: when I test for the significant contribution of monomials (forward elimination) anova(rda(Helling ~ I(x^2)+Condition(x)+Condition(y))) performs the permutation test as expected, whereas anova(rda(Helling ~ I(y^2)+Condition(x)+Condition(y))) Returns this error message: Error in names-.default(`*tmp*`, value = Model) : attempt to set an attribute on NULL Thanks again for your help Kind regards, Helene Helene MORLON University of California, Merced -Original Message- From: Jari Oksanen [mailto:[EMAIL PROTECTED] Sent: Monday, February 26, 2007 11:27 PM To: r-help@stat.math.ethz.ch Cc: [EMAIL PROTECTED] Subject: [R] RDA and trend surface regression 'm performing RDA on plant presence/absence data, constrained by geographical locations. I'd like to constrain the RDA by the extended matrix of geographical coordinates -ie the matrix of geographical coordinates completed by adding all terms of a cubic trend surface regression- . This is the command I use (package vegan): rda(Helling ~ x+y+x*y+x^2+y^2+x*y^2+y*x^2+x^3+y^3) where Helling is the matrix of Hellinger-transformed presence/absence data The result returned by R is exactly the same as the one given by: anova(rda(Helling ~ x+y) Ie the quadratic and cubic terms are not taken into account You must *I*solate the polynomial terms with function I (AsIs) so that they are not interpreted as formula operators: rda(Helling ~ x + y + I(x*y) + I(x^2) + I(y^2) + I(x*y^2) + I(y*x^2) + I(x^3) + I(y^3)) If you don't have the interaction terms, then it is easier and better (numerically) to use poly(): rda(Helling ~ poly(x, 3) + poly(y, 3)) Another issue is that in my opinion using polynomial constraints is an Extremely Bad Idea(TM). cheers, Jari Oksanen __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] RDA and trend surface regression
Helene, My point was only that RDA may fit a quadratic model for the terms specified in your model. The terms that you had specified were already higher order polynomials (some cubic). So a QDA classifier with the model terms that you specified my be a fifth order polynomial in the original data. I don't know the reference you cite or even the subject-matter specifics. I'm just a simple cave man (for you SNL fans). But I do know that there are more reliable ways to get nonlinear classification boundaries than using x^5. If you want a quadratic model, I would suggest that you use QDA with the predictors in the original units (or see Hastie's book for a good example of using higher order terms with LDA). Looking at your email, you want a a variation partitioning analyses. RDA works best as a classification technique. Perhaps a multivariate ANOVA model may be a more direct way to meet your needs. There is a connection between LDA and some multivariate linear models, but I don't know of a similar connection to RDA. Max -Original Message- From: MORLON [mailto:[EMAIL PROTECTED] Sent: Tuesday, February 27, 2007 12:53 PM To: 'Jari Oksanen'; r-help@stat.math.ethz.ch Cc: Kuhn, Max Subject: RE: [R] RDA and trend surface regression Thanks a lot for your answers, I am concerned by your advice not to use polynomial constraints, or to use QDA instead of RDA. My final goal is to perform variation partitioning using partial RDA to assess the relative importance of environmental vs spatial variables. For the spatial analyses, trend surface analysis (polynomial constraints) is recommended in Legendre and Legendre 1998 (p739). Is there a better method to integrate space as an explanatory variable in a variation partitioning analyses? Also, I don't understand this: when I test for the significant contribution of monomials (forward elimination) anova(rda(Helling ~ I(x^2)+Condition(x)+Condition(y))) performs the permutation test as expected, whereas anova(rda(Helling ~ I(y^2)+Condition(x)+Condition(y))) Returns this error message: Error in names-.default(`*tmp*`, value = Model) : attempt to set an attribute on NULL Thanks again for your help Kind regards, Helene Helene MORLON University of California, Merced -Original Message- From: Jari Oksanen [mailto:[EMAIL PROTECTED] Sent: Monday, February 26, 2007 11:27 PM To: r-help@stat.math.ethz.ch Cc: [EMAIL PROTECTED] Subject: [R] RDA and trend surface regression 'm performing RDA on plant presence/absence data, constrained by geographical locations. I'd like to constrain the RDA by the extended matrix of geographical coordinates -ie the matrix of geographical coordinates completed by adding all terms of a cubic trend surface regression- . This is the command I use (package vegan): rda(Helling ~ x+y+x*y+x^2+y^2+x*y^2+y*x^2+x^3+y^3) where Helling is the matrix of Hellinger-transformed presence/absence data The result returned by R is exactly the same as the one given by: anova(rda(Helling ~ x+y) Ie the quadratic and cubic terms are not taken into account You must *I*solate the polynomial terms with function I (AsIs) so that they are not interpreted as formula operators: rda(Helling ~ x + y + I(x*y) + I(x^2) + I(y^2) + I(x*y^2) + I(y*x^2) + I(x^3) + I(y^3)) If you don't have the interaction terms, then it is easier and better (numerically) to use poly(): rda(Helling ~ poly(x, 3) + poly(y, 3)) Another issue is that in my opinion using polynomial constraints is an Extremely Bad Idea(TM). cheers, Jari Oksanen -- LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}} __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] RDA and trend surface regression
On Tue, 2007-02-27 at 13:13 -0500, Kuhn, Max wrote: Helene, My point was only that RDA may fit a quadratic model for the terms specified in your model. The terms that you had specified were already higher order polynomials (some cubic). So a QDA classifier with the model terms that you specified my be a fifth order polynomial in the original data. I don't know the reference you cite or even the subject-matter specifics. I'm just a simple cave man (for you SNL fans). But I do know that there are more reliable ways to get nonlinear classification boundaries than using x^5. I doubt that Helene is trying to do a classification - unless you consider classification to mean that all rows/samples are in different groups (i.e. n samples therefore n groups) - which is how RDA (Redundancy Analysis) is used in ecology. You could take a look at multispati in package ade4 for a different way to handle spatial constraints. There is also the principle coordinates analysis of neighbour matrices (PCNM) method - not sure this is coded anywhere in R yet though. Here are two references that may be useful: Dray, S., P. Legendre, and P. R. Peres- Neto. 2006. Spatial modeling: a comprehensive framework for principal coordinate analysis of neighbor matrices (PCNM). Ecological Modelling, in press. Griffith, D. A., and P. R. Peres- Neto. 2006. Spatial modeling in ecology: the flexibility of eigenfunction spatial analyses. Ecology, in press. HTH G If you want a quadratic model, I would suggest that you use QDA with the predictors in the original units (or see Hastie's book for a good example of using higher order terms with LDA). Looking at your email, you want a a variation partitioning analyses. RDA works best as a classification technique. Perhaps a multivariate ANOVA model may be a more direct way to meet your needs. There is a connection between LDA and some multivariate linear models, but I don't know of a similar connection to RDA. Max -Original Message- From: MORLON [mailto:[EMAIL PROTECTED] Sent: Tuesday, February 27, 2007 12:53 PM To: 'Jari Oksanen'; r-help@stat.math.ethz.ch Cc: Kuhn, Max Subject: RE: [R] RDA and trend surface regression Thanks a lot for your answers, I am concerned by your advice not to use polynomial constraints, or to use QDA instead of RDA. My final goal is to perform variation partitioning using partial RDA to assess the relative importance of environmental vs spatial variables. For the spatial analyses, trend surface analysis (polynomial constraints) is recommended in Legendre and Legendre 1998 (p739). Is there a better method to integrate space as an explanatory variable in a variation partitioning analyses? Also, I don't understand this: when I test for the significant contribution of monomials (forward elimination) anova(rda(Helling ~ I(x^2)+Condition(x)+Condition(y))) performs the permutation test as expected, whereas anova(rda(Helling ~ I(y^2)+Condition(x)+Condition(y))) Returns this error message: Error in names-.default(`*tmp*`, value = Model) : attempt to set an attribute on NULL Thanks again for your help Kind regards, Helene Helene MORLON University of California, Merced -Original Message- From: Jari Oksanen [mailto:[EMAIL PROTECTED] Sent: Monday, February 26, 2007 11:27 PM To: r-help@stat.math.ethz.ch Cc: [EMAIL PROTECTED] Subject: [R] RDA and trend surface regression 'm performing RDA on plant presence/absence data, constrained by geographical locations. I'd like to constrain the RDA by the extended matrix of geographical coordinates -ie the matrix of geographical coordinates completed by adding all terms of a cubic trend surface regression- . This is the command I use (package vegan): rda(Helling ~ x+y+x*y+x^2+y^2+x*y^2+y*x^2+x^3+y^3) where Helling is the matrix of Hellinger-transformed presence/absence data The result returned by R is exactly the same as the one given by: anova(rda(Helling ~ x+y) Ie the quadratic and cubic terms are not taken into account You must *I*solate the polynomial terms with function I (AsIs) so that they are not interpreted as formula operators: rda(Helling ~ x + y + I(x*y) + I(x^2) + I(y^2) + I(x*y^2) + I(y*x^2) + I(x^3) + I(y^3)) If you don't have the interaction terms, then it is easier and better (numerically) to use poly(): rda(Helling ~ poly(x, 3) + poly(y, 3)) Another issue is that in my opinion using polynomial constraints is an Extremely Bad Idea(TM). cheers, Jari Oksanen -- LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}} __ 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] RDA and trend surface regression
On Tue, 27 Feb 2007, Gavin Simpson wrote: On Tue, 2007-02-27 at 13:13 -0500, Kuhn, Max wrote: Helene, My point was only that RDA may fit a quadratic model for the terms specified in your model. The terms that you had specified were already higher order polynomials (some cubic). So a QDA classifier with the model terms that you specified my be a fifth order polynomial in the original data. I don't know the reference you cite or even the subject-matter specifics. I'm just a simple cave man (for you SNL fans). But I do know that there are more reliable ways to get nonlinear classification boundaries than using x^5. I doubt that Helene is trying to do a classification - unless you consider classification to mean that all rows/samples are in different groups (i.e. n samples therefore n groups) - which is how RDA (Redundancy Analysis) is used in ecology. You could take a look at multispati in package ade4 for a different way to handle spatial constraints. There is also the principle coordinates analysis of neighbour matrices (PCNM) method - not sure this is coded anywhere in R yet though. Here are two references that may be useful: Dray, S., P. Legendre, and P. R. Peres- Neto. 2006. Spatial modeling: a comprehensive framework for principal coordinate analysis of neighbor matrices (PCNM). Ecological Modelling, in press. Griffith, D. A., and P. R. Peres- Neto. 2006. Spatial modeling in ecology: the flexibility of eigenfunction spatial analyses. Ecology, in press. Pedro Peres-Neto helped cast his matlab original to R as function ME() in the spdep package, at least partly to see if it worked like SpatialFiltering() in spdep, based on a forthcoming paper by Tiefelsdorf and Griffith; code suggestions from Stéphane Dray were also used. For sample data sets, including those used in the papers, ME() reproduces the original results, and reproduces results from the matlab code from which it was derived with possible differences for the stopping rule of the semi-parametric stage - the Oribatid mites data set is in ade4. Roger HTH G If you want a quadratic model, I would suggest that you use QDA with the predictors in the original units (or see Hastie's book for a good example of using higher order terms with LDA). Looking at your email, you want a a variation partitioning analyses. RDA works best as a classification technique. Perhaps a multivariate ANOVA model may be a more direct way to meet your needs. There is a connection between LDA and some multivariate linear models, but I don't know of a similar connection to RDA. Max -Original Message- From: MORLON [mailto:[EMAIL PROTECTED] Sent: Tuesday, February 27, 2007 12:53 PM To: 'Jari Oksanen'; r-help@stat.math.ethz.ch Cc: Kuhn, Max Subject: RE: [R] RDA and trend surface regression Thanks a lot for your answers, I am concerned by your advice not to use polynomial constraints, or to use QDA instead of RDA. My final goal is to perform variation partitioning using partial RDA to assess the relative importance of environmental vs spatial variables. For the spatial analyses, trend surface analysis (polynomial constraints) is recommended in Legendre and Legendre 1998 (p739). Is there a better method to integrate space as an explanatory variable in a variation partitioning analyses? Also, I don't understand this: when I test for the significant contribution of monomials (forward elimination) anova(rda(Helling ~ I(x^2)+Condition(x)+Condition(y))) performs the permutation test as expected, whereas anova(rda(Helling ~ I(y^2)+Condition(x)+Condition(y))) Returns this error message: Error in names-.default(`*tmp*`, value = Model) : attempt to set an attribute on NULL Thanks again for your help Kind regards, Helene Helene MORLON University of California, Merced -Original Message- From: Jari Oksanen [mailto:[EMAIL PROTECTED] Sent: Monday, February 26, 2007 11:27 PM To: r-help@stat.math.ethz.ch Cc: [EMAIL PROTECTED] Subject: [R] RDA and trend surface regression 'm performing RDA on plant presence/absence data, constrained by geographical locations. I'd like to constrain the RDA by the extended matrix of geographical coordinates -ie the matrix of geographical coordinates completed by adding all terms of a cubic trend surface regression- . This is the command I use (package vegan): rda(Helling ~ x+y+x*y+x^2+y^2+x*y^2+y*x^2+x^3+y^3) where Helling is the matrix of Hellinger-transformed presence/absence data The result returned by R is exactly the same as the one given by: anova(rda(Helling ~ x+y) Ie the quadratic and cubic terms are not taken into account You must *I*solate the polynomial terms with function I (AsIs) so that they are not interpreted
Re: [R] RDA and trend surface regression
On 27 Feb 2007, at 20:55, Gavin Simpson wrote: On Tue, 2007-02-27 at 13:13 -0500, Kuhn, Max wrote: Helene, My point was only that RDA may fit a quadratic model for the terms specified in your model. The terms that you had specified were already higher order polynomials (some cubic). So a QDA classifier with the model terms that you specified my be a fifth order polynomial in the original data. I don't know the reference you cite or even the subject-matter specifics. I'm just a simple cave man (for you SNL fans). But I do know that there are more reliable ways to get nonlinear classification boundaries than using x^5. I doubt that Helene is trying to do a classification - unless you consider classification to mean that all rows/samples are in different groups (i.e. n samples therefore n groups) - which is how RDA (Redundancy Analysis) is used in ecology. You could take a look at multispati in package ade4 for a different way to handle spatial constraints. There is also the principle coordinates analysis of neighbour matrices (PCNM) method - not sure this is coded anywhere in R yet though. Here are two references that may be useful: Stéphane Dray has R code for finding PCNM matrices. Google for his name: it's not that common. I also have a copy of his function and can send it if really needed, but it may be better to check Dray's page first. Stéphane Dray says think that not all functions need be in CRAN. May be true, but I think it might help many people. There are at least three reasons why not use polynomial constraints in RDA. Max Kuhn mentioned one: polynomials typically flip wildly at margins (or they are unstable in more neutral speech). Second reason is that they are almost impossible to interpret in ordination display. The third reason is that RDA (or CCA) avoid some ordination artefacts (curvature, horseshoe, arc etc.) just because the constraints are linear: allowing them to be curved allows curved solutions. These arguments are not necessarily valid if you only want to have variance partitioning, or if you use polynomial conditions (partial out polynomial effects in Canoco language). In that case it may make sense to use quadratic (or polynomial) constraints or conditions. cheers, Jari Oksanen __ 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 and provide commented, minimal, self-contained, reproducible code.
[R] RDA and trend surface regression
Dear all, I'm performing RDA on plant presence/absence data, constrained by geographical locations. I'd like to constrain the RDA by the extended matrix of geographical coordinates -ie the matrix of geographical coordinates completed by adding all terms of a cubic trend surface regression- . This is the command I use (package vegan): rda(Helling ~ x+y+x*y+x^2+y^2+x*y^2+y*x^2+x^3+y^3) where Helling is the matrix of Hellinger-transformed presence/absence data The result returned by R is exactly the same as the one given by: anova(rda(Helling ~ x+y) Ie the quadratic and cubic terms are not taken into account I hope you can help me with that: how can I perform a RDA on an extended matrix of geographical coordinates in R?. Thank you very much in advance, Helene Morlon University of California, Merced [EMAIL PROTECTED] [[alternative HTML version deleted]] __ 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 and provide commented, minimal, self-contained, reproducible code.
Re: [R] RDA and trend surface regression
Helene, You will have to give us more information, such as your system/versions and a small reproducible example. We try to stress that questions are more easily answered when there are a lot of specific details given and a reproducible case can be tested. Here are two comments though: 1. The quadratic terms probably are not showing up because you are not using a proper model formula for the task. See: http://cran.r-project.org/doc/manuals/R-intro.html#Formulae-for-statisti cal-models Specifically, the part that says I(M): Insulate M. Inside M all operators have their normal arithmetic meaning, and that term appears in the model matrix. is important. So, as an example from ?rda: x - rda(Species ~ (Sepal.Length+Sepal.Width)^2 + Sepal.Width^2, data = iris) would not work for the squared term, but x - rda(Species ~ (Sepal.Length+Sepal.Width)^2 + I(Sepal.Width^2), data = iris) would. 2. RDA is fitting models at or between LDA and QDA. So a QDA model with quadratic terms would be quartic discriminant analysis. Of course, there are no rules against this, but high order polynomials can do weird things in the tail (which would be the edges of the space defined by your training data). If your data are that nonlinear, there are much better ways of classifying data. I'd suggests getting a copy of Hastie et all (2001) or MASS. Max -Original Message- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of MORLON Sent: Monday, February 26, 2007 7:14 PM To: r-help@stat.math.ethz.ch Subject: [R] RDA and trend surface regression Dear all, I'm performing RDA on plant presence/absence data, constrained by geographical locations. I'd like to constrain the RDA by the extended matrix of geographical coordinates -ie the matrix of geographical coordinates completed by adding all terms of a cubic trend surface regression- . This is the command I use (package vegan): rda(Helling ~ x+y+x*y+x^2+y^2+x*y^2+y*x^2+x^3+y^3) where Helling is the matrix of Hellinger-transformed presence/absence data The result returned by R is exactly the same as the one given by: anova(rda(Helling ~ x+y) Ie the quadratic and cubic terms are not taken into account I hope you can help me with that: how can I perform a RDA on an extended matrix of geographical coordinates in R?. Thank you very much in advance, Helene Morlon University of California, Merced [EMAIL PROTECTED] [[alternative HTML version deleted]] __ 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 and provide commented, minimal, self-contained, reproducible code. -- LEGAL NOTICE\ Unless expressly stated otherwise, this messag...{{dropped}} __ 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 and provide commented, minimal, self-contained, reproducible code.
[R] RDA and trend surface regression
'm performing RDA on plant presence/absence data, constrained by geographical locations. I'd like to constrain the RDA by the extended matrix of geographical coordinates -ie the matrix of geographical coordinates completed by adding all terms of a cubic trend surface regression- . This is the command I use (package vegan): rda(Helling ~ x+y+x*y+x^2+y^2+x*y^2+y*x^2+x^3+y^3) where Helling is the matrix of Hellinger-transformed presence/absence data The result returned by R is exactly the same as the one given by: anova(rda(Helling ~ x+y) Ie the quadratic and cubic terms are not taken into account You must *I*solate the polynomial terms with function I (AsIs) so that they are not interpreted as formula operators: rda(Helling ~ x + y + I(x*y) + I(x^2) + I(y^2) + I(x*y^2) + I(y*x^2) + I(x^3) + I(y^3)) If you don't have the interaction terms, then it is easier and better (numerically) to use poly(): rda(Helling ~ poly(x, 3) + poly(y, 3)) Another issue is that in my opinion using polynomial constraints is an Extremely Bad Idea(TM). cheers, Jari Oksanen __ 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 and provide commented, minimal, self-contained, reproducible code.