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.
