I am not sure that method works. It appears to be doing something close, but returns too many slope coefficients, since I think it is returning interaction terms of degree smaller and greater than what was passed to it.
Here is a small example of degree 2: > X = data.frame(cbind(c(1,2,3), c(4,5,6))) > X X1 X2 1 1 4 2 2 5 3 3 6 > y = c(0,1,0) > mymodel2 = glm(y ~ poly(X$X1,2)*poly(X$X2,2), > family=binomial(link="logit")) # We see slope coefficients returned for interaction terms of degree 1, 3, and 4. > summary(mymodel2) Estimate Std. Error z value Pr(>|z|) (Intercept) -7.855e+00 4.588e+04 0 1 poly(X$X1, 2)1 6.059e-15 7.946e+04 0 1 poly(X$X1, 2)2 -3.848e+01 7.946e+04 0 1 poly(X$X2, 2)1 NA NA NA NA poly(X$X2, 2)2 NA NA NA NA poly(X$X1, 2)1:poly(X$X2, 2)1 NA NA NA NA poly(X$X1, 2)2:poly(X$X2, 2)1 NA NA NA NA poly(X$X1, 2)1:poly(X$X2, 2)2 NA NA NA NA poly(X$X1, 2)2:poly(X$X2, 2)2 NA NA NA NA # Data used in the model > mymodel2$model y poly(X$X1, 2).1 poly(X$X1, 2).2 poly(X$X2, 2).1 poly(X$X2, 2).2 1 0 -7.071068e-01 4.082483e-01 -7.071068e-01 4.082483e-01 2 1 4.350720e-18 -8.164966e-01 4.350720e-18 -8.164966e-01 3 0 7.071068e-01 4.082483e-01 7.071068e-01 4.082483e-01 where instead I would expect the data to be like 1 4 16 4 10 25 9 18 36 Any other ideas? Many thanks. -- View this message in context: http://r.789695.n4.nabble.com/Specify-model-with-polynomial-interaction-terms-up-to-degree-n-tp4635130p4635196.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help@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.