Currently, if any data vectors given to poly() contain missing values then poly throws an error. There is no need for this if the coefs argument is supplied nor if raw=TRUE - it could just return a row of NA's for the elements of the data containing NA's. The current behavior causes unneeded errors in predict when newdata contains missing values: > lmFit <- lm(data=data.frame(x=1:7, y=sin(1:7)), y ~ poly(x, 3)) > predict(lmFit, newdata=list(x=c(1,2,3))) 1 2 3 0.9324742 0.7177399 0.1104601 > predict(lmFit, newdata=list(x=c(1,2,3,NA))) Error in poly(x, 3, coefs = list(alpha = c(4, 4, 4), norm2 = c(1, 7, 28, : missing values are not allowed in 'poly' I think it would be better for the last to return 1 2 3 4 0.9324742 0.7177399 0.1104601 NA Checking for NA's in the orthogonal polynomial branch of the function only when is.null(coefs) would fix the problem. It is hard to enough to get people to use poly(x,3) instead of I(x)+I(x^2)+I(x^3) without this problem.

You are entirely right. I'm committing a patch for this to both R-devel and R 3.2.2 patched which will become R 3.2.3 in about one month. Thank you very much, Bill!