[R] Specify model with polynomial interaction terms up to degree n

[YTP]

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.