Solved – How would you report (in publication) the results of a linear model fit using the poly function in R

polynomialrregressionregression coefficients

@John recently pointed out to me that R's poly function produces less correlated values (more orthogonal) to fit polynomial predictors, i.e. the transformed predictors have a lower correlation with each other than if I just mean centered the 1st order predictor prior to generating the nth order predictor(s). R's poly accomplishes this using an algorithm from Kennedy, W. J. Jr and Gentle, J. E. (1980) Statistical Computing Marcel Dekker, pg 343-344. I checked the book out from the library, but I have to say that I am not quite following what is going on. I've also made some plots comparing these values to the raw values and I sort of can see what is going on, but I still don't I fully understand it. It seems like the data is being processed and subjected to some transformation that involves the mean of the values and some prediction coefficients. I can always use predict.poly to turn the values from real values to poly values. My question is how should I report the results of a model fit with poly? The raw slope coefficients refer to a transformed value, if I report them, then I assume I should also report some of the coefficients from poly itself as well. Which ones? In any particular format? Is anybody going to really understand them? Is this a compelling reason to stick to raw mean-centered coefficients?

Best Answer

With modern matrix algebra software, I feel that orthogonal polynomials get in the way more than they help. I prefer to use ordinary polynomials for this reason, for example the pol function in the R rms package.

Related Solutions

Solved – get wildly different results for poly(raw=T) vs. poly()

I think this is a bug in the predict function (and hence my fault), which in fact nlme does not share. (Edit: should be fixed in most recent R-forge version of lme4.) See below for an example ...

I think your understanding of orthogonal polynomials is probably just fine. The tricky thing you need to know about them if you are trying to write a predict method for a class of models is that the basis for the orthogonal polynomials is defined based on a given set of data, so if you naively (like I did!) use model.matrix to try to generate the design matrix for a new set of data, you get a new basis -- which no longer makes sense with the old parameters. Until I get this fixed, I may need to put a trap in that tells people that predict doesn't work with orthogonal polynomial bases (or spline bases, which have the same property).

d <- expand.grid(x=seq(0,1,length=50),f=LETTERS[1:10])
set.seed(1001)
u.int <- rnorm(10,sd=0.5)
u.slope <- rnorm(10,sd=0.2)
u.quad <- rnorm(10,sd=0.1)
d <- transform(d,
               ypred = (1+u.int[f])+
               (2+u.slope[f])*x-
               (1+u.quad[f])*x^2)
d$y <- rnorm(nrow(d),mean=d$ypred,sd=0.2)
ggplot(d,aes(x=x,y=y,colour=f))+geom_line()+
    geom_line(aes(y=ypred),linetype=2)

library(lme4)
fm1 <- lmer(y~poly(x,2,raw=TRUE)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)


fm2 <- lmer(y~poly(x,2)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)
newdat <- data.frame(x=unique(d$x))
plot(predict(fm1,newdata=newdat,REform=NA))
lines(predict(fm2,newdata=newdat,REform=NA),col=2)
detach("package:lme4")

library(nlme)
fm3 <- lme(y~poly(x,2,raw=TRUE),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)
VarCorr(fm3)

fm4 <- lme(y~poly(x,2),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)

newdat <- data.frame(x=unique(d$x))
lines(predict(fm3,newdata=newdat,level=0),col=4)
lines(predict(fm4,newdata=newdat,level=0),col=5)

Solved – How to present results of a Lasso using glmnet

My understanding is that you can't necessarily say much about which variables are "important" or have "real" effects based on whether their coefficients are nonzero. To give an extreme example, if you have two predictors that are perfectly collinear, the lasso will pick one of them essentially at random to get the full weight and the other one will get zero weight.

This paper, which includes one of the authors of glmnet, presents some glmnet-based analyses (see especially: the Introduction, Sections 2.3 and 4.3, and Tables 4 and 5). Glancing through, it looks like they didn't calculate P-valued directly from the glmnet model. They did calculate two different kinds of P-values using other methods, but it doesn't look like they fully trust either of them.

I'm not 100% sure what you're suggesting in terms of plotting methods, but I think it sounds reasonable.

Hope that helps.

Related Question