4
$\begingroup$

Note: I've updated the example case code, there were some errors in the previous version

Cross posted to R-help, because I half suspect this is 'unexpected behaviour'.

I want to predict values from an existing lm (linear model, e.g. lm.obj) result in R using a new set of predictor variables (e.g. newdata). Specifically, I am interested in the predicted y value at the mean, 1 SD above of the mean, and 1 SD below the mean for each predictor. However, it seems that because my linear models was made by calling scale() on the target predictor that predict exits with an error, "Error in scale(xxA, center = 9.7846094491829, scale = 0.959413568556403) : object 'xxA' not found". By debugging predict, I can see that the error occurs in a call to model.frame. By debugging model frame I can see the error occurs with this command: variables <- eval(predvars, data, env); it seems likely that the error is because predvars looks like this:

list(scale(xxA, center = 10.2058714830537, scale = 0.984627257169526), scale(xxB, center = 20.4491690881149, scale = 1.13765718273923))

An example case:

   dat <- data.frame(xxA = rnorm(20,10), xxB = rnorm(20,20))
   dat$out <- with(dat,xxA+xxB+xxA*xxB+rnorm(20,20))
   lm.res.scale <- lm(out ~ scale(xxA)*scale(xxB),data=dat)
   my.data <- lm.res.scale$model #load the data from the lm object
   newdata <- expand.grid(X1=c(-1,0,1),X2=c(-1,0,1))
   names(newdata) <- c("scale(xxA)","scale(xxB)")
   newdata$Y <- predict(lm.res.scale,newdata)

Is there something I could do before passing newdata or lm.obj to predict() that would prevent the error? I tried:

From the help file it looks like I might be able to do something with the terms, argument but I haven't quite figured out what I would need to do. Alternatively, is there a fix for model.frame that would prevent the error? Should predict() behave this way?

Additional Details: However, I really want a solution that, in one step will provide values like:

coef(lm.res.scale)[1]+
coef(lm.res.scale)[2]*newdata[,1]+
coef(lm.res.scale)[3]*newdata[,2]+
coef(lm.res.scale)[4]*newdata[,1]*newdata[,2]

I think that should be exactly what predict() should do. That is, I think my example code should be equivalent to:

   dat <- data.frame(xxA = rnorm(20,10), xxB = rnorm(20,20))
   dat$out <- with(dat,xxA+xxB+xxA*xxB+rnorm(20,20))
   #rescaling outside of lm
   X1 <- with(dat,as.vector(scale(xxA)))
   X2 <- with(dat,as.vector(scale(xxB)))
   y  <- with(dat,out)
   lm.res.correct <- lm(y~X1*X2)
   my.data <- lm.res.correct$model #load the data from the lm object
   newdata <- expand.grid(X1=c(-1,0,1),X2=c(-1,0,1))
   #No need to rename newdata as it matches my lm object already
   newdata$Y <- predict(lm.res.correct,newdata)

Notably, adjusting my formula to include as.vector() does not solve the problem with my attempt to use predict() directly with newdata.

$\endgroup$
  • $\begingroup$ these types of questions might get better answers in stackoverflow.com. $\endgroup$ – mpiktas Jan 29 '11 at 8:16
  • $\begingroup$ If xxA and xxB are standardized before calling lm(), it seems to work. E.g., dat <- cbind.data.frame(dat, xxAstd=scale(dat$xxA), xxBstd=scale(dat$xxB)); lm.res.scale <- lm(out ~ xxAstd*xxBstd,data=dat); names(newdata) <- c("xxAstd","xxBstd"); newdata$Y <- predict(lm.res.scale,newdata). Not sure it's what you seek, though. $\endgroup$ – chl Jan 29 '11 at 11:50
  • $\begingroup$ Yeah chl, I noticed that, but was disappointed I couldn't call scale from within lm and have predict work. $\endgroup$ – russellpierce Jan 29 '11 at 16:34
5
$\begingroup$

When you use the predict with newdata argument you must supply the the data.frame with the same column names. In your code you have

newdata <- expand.grid(X1=c(-1,0,1),X2=c(-1,0,1))
names(newdata) <- c("scale(xxA)","scale(xxB)")

But the formula supplied to lm object is

out ~ scale(xxA)*scale(xxB)

So when you call the predict, it tries to find objects xxA and xxB in your data and apply function scale as per your initial request. But all R finds are objects scale(xxA) and scale(xxB). So naturally it produces the error.

Now if you supply correctly named newdata

newdata <- expand.grid(xxA=c(-1,0,1),xxB=c(-1,0,1)) 

and try to use it for prediction

newdata$Y <- predict(lm.res.scale,newdata)

R will remember the how it scaled original data and apply the same scaling to your new data. In this case supplied value -1 for xxA will be subtracted the original mean of xxA and divided by the original standard value of xxA. If you want to get prediction of 1 S.D below the mean, you will need to supply this value. In your case then newdata should look like this:

newdata <- expand.grid(xxA=mean(dat$xxA)+sd(dat$xxA)*c(-1,0,1),xxB=mean(dat$xxB)+sd(dat$xxB)*c(-1,0,1))

I gathered all the solutions in one place to compare:

##Prepare data
set.seed(1)
dat <- data.frame(xxA = rnorm(20,10), xxB = rnorm(20,20))
dat$out <- with(dat,xxA+xxB+xxA*xxB+rnorm(20,20))
dat <- within(dat,{
    X1 <- as.numeric(scale(xxA))
    X2 <- as.numeric(scale(xxB))
    })

##Estimate the models
lm.res.scale <- lm(out ~ scale(xxA)*scale(xxB),data=dat)
lm.res.correct <- lm(out~X1*X2,data=dat)
lm.mod <- lm(out ~ I(scale(xxA))*I(scale(xxB)), data=dat)
rms.res <- Glm(out ~ scale(xxA)*scale(xxB),data=dat)


##Build data for prediction
newdata <- expand.grid(xxA=c(-1,0,1),xxB=c(-1,0,1))
newdata$X1<-newdata$xxA
newdata$X2<-newdata$xxB

##Gather the predictions
newdata$Yscaled <- predict(lm.res.scale,newdata)
newdata$Ycorrect <- predict(lm.res.correct,newdata)
newdata$YwithI <- predict(lm.mod,newdata)
newdata$Ywithrms <- Predict(rms.res,xxA=c(-1,0,1),xxB=c(-1,0,1),conf.int=FALSE)[,3]

##Build alternative data for prediction
newdata2 <- expand.grid(xxA=mean(dat$xxA)+sd(dat$xxA)*c(-1,0,1),xxB=mean(dat$xxB)+sd(dat$xxB)*c(-1,0,1))

#Predict
newdata$Yorigsc <- predict(lm.res.scale,newdata2)

I used set.seed so the results should be the same if you try to repeat it. The newdata looks like this:

> newdata
  xxA xxB X1 X2  Yscaled Ycorrect   YwithI Ywithrms  Yorigsc
1  -1  -1 -1 -1 25.79709 225.9562 221.7517 221.7517 225.9562
2   0  -1  0 -1 25.63030 244.5181 243.0404 243.0404 244.5181
3   1  -1  1 -1 25.46351 263.0800 264.3291 264.3291 263.0800
4  -1   0 -1  0 25.36341 234.6981 231.7012 231.7012 234.6981
5   0   0  0  0 26.21499 254.0704 254.0704 254.0704 254.0704
6   1   0  1  0 27.06657 273.4427 276.4396 276.4396 273.4427
7  -1   1 -1  1 24.92972 243.4400 241.6507 241.6507 243.4400
8   0   1  0  1 26.79967 263.6227 265.1004 265.1004 263.6227
9   1   1  1  1 28.66962 283.8054 288.5501 288.5501 283.8054

As expected Yscaled produces not the result we need since the original scaling is applied. In the case when we scale data before lm (Ycorrect) and when we supply alternative unscaled values (Yorigsc) results coincide and are the ones needed.

Now the other prediction methods give different results. This happens since R is forced to forget the original scaling using formula

out ~ I(scale(xxA))*I(scale(xxB))

or package rms. But when we use predict, the values are still scaled, but now according to supplied values of xxA and xxB. This is best illustrated by following statement, which in some way mimics what predict does with the data:

> eval(expression(cbind(scale(xxA),scale(xxB))),env=as.list(newdata))
           [,1]      [,2]
 [1,] -1.154701 -1.154701
 [2,]  0.000000 -1.154701
 [3,]  1.154701 -1.154701
 [4,] -1.154701  0.000000
 [5,]  0.000000  0.000000
 [6,]  1.154701  0.000000
 [7,] -1.154701  1.154701
 [8,]  0.000000  1.154701
 [9,]  1.154701  1.154701

We can see that in this case, scaling does not change original values too much, but this is even worse, since the values from predict look reasonable, when in fact they are wrong.

$\endgroup$
  • $\begingroup$ mean(newdata[,"Y"]); mean(dat[,"out:). If predict was working as expected these values should be near one another. They aren't. $\endgroup$ – russellpierce Jan 30 '11 at 16:24
  • $\begingroup$ What does, "this happens since R is forced to forget the original scaling using formula", mean? What is meant by saying that R knows or does not know the original scaling? $\endgroup$ – russellpierce Jan 31 '11 at 1:05
  • $\begingroup$ Ah... I see - you mean that the predictors in new data themselves are subjected to Z scaling prior to the generation of the prediction? $\endgroup$ – russellpierce Jan 31 '11 at 1:29
  • $\begingroup$ ... so I think the theme here that there is no way to use the scale command within the call to lm and pass predict the Z values I want and have everything turn out smoothy in one step the way I might hope. $\endgroup$ – russellpierce Jan 31 '11 at 1:37
1
$\begingroup$

See if this console scrape shows what you might have wanted. Changed the names of the new data columns:

> dat <- data.frame(xxA = rnorm(20,10), xxB = rnorm(10,20))
>    dat$out <- with(dat,xxA+xxB+xxA*xxB+rnorm(20,20))
>    xVar <- "scale(xxA)"
>    traceVar <- "scale(xxB)"
>    DVname <- "out"
>    lm.obj <- lm.res.scale <- lm(out ~ scale(xxA)*scale(xxB),data=dat)
>    my.data <- lm.obj$model #load the data from the lm object
>    X1 <- my.data[,xVar]
>    X2 <- my.data[,traceVar]
>    DV <- lm.obj$model[,DVname]
>    newdata <- expand.grid(X1=c(-1,0,1),X2=c(-1,0,1))
>    newdata$xxA <- newdata$X1 * sd(my.data[,xVar])
>    newdata$xxB <- newdata$X2 * sd(my.data[,traceVar])
>    names(newdata) <- names(dat) #have to rename to original variable names for predict to work
>    newdata$Y <- predict(lm.obj,newdata)
> str(newdata)
'data.frame':   9 obs. of  5 variables:
 $ xxA: num  -1 0 1 -1 0 1 -1 0 1
     $ xxB: num  -1 -1 -1 0 0 0 1 1 1
 $ out: num  -1 0 1 -1 0 1 -1 0 1
     $ NA : num  -1 -1 -1 0 0 0 1 1 1
 $ Y  : num  47.7 44.7 41.8 46.5 44.7 ...
     - attr(*, "out.attrs")=List of 2
      ..$ dim     : Named int  3 3
  .. ..- attr(*, "names")= chr  "X1" "X2"
  ..$ dimnames:List of 2
  .. ..$ X1: chr  "X1=-1" "X1= 0" "X1= 1"
  .. ..$ X2: chr  "X2=-1" "X2= 0" "X2= 1"

EDIT With your requested revisions, this is what I would suggest:

newdat <- expand.grid(X1=c(-1,0,1),X2=c(-1,0,1))
newdat$xxA <- mean(dat$xxA) + newdat$X1*sd(dat$xxA)
newdat$xxB <- mean(dat$xxB) + newdat$X2*sd(dat$xxB)
 newdata$Y <- predict(lm.res.scale,newdat[, c("xxA","xxB") ] )
 newdat$Y <- predict(lm.res.scale,newdat[, c("xxA","xxB") ] )
 newdat

  X1 X2       xxA      xxB        Y
1 -1 -1  8.931712 18.59229 214.0096
2  0 -1  9.830117 18.59229 231.4839
3  1 -1 10.728522 18.59229 248.9583
4 -1  0  8.931712 19.77820 225.3516
5  0  0  9.830117 19.77820 243.9916
6  1  0 10.728522 19.77820 262.6315
7 -1  1  8.931712 20.96410 236.6936
8  0  1  9.830117 20.96410 256.4992
9  1  1 10.728522 20.96410 276.3047

Notice you don't need to do all those back-flips to match the model terms inside the lm-object. You only need to provide it with properly named data on the scale of the original data. The predict machinery will take care of the scaling.

$\endgroup$
  • $\begingroup$ If I use names(dat) instead of c(xVar,traceVar) when renaming newdata prior to using predict, then predict doesn't do what I'd expect. I'm sorry that my previous example code was unclear. Now, I've provided a more clear example above. What I want is to find the predicted y value at the mean, 1 SD above of the mean, and 1 SD below the mean for each predictor. $\endgroup$ – russellpierce Jan 29 '11 at 1:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.