R help - Apr 2008 - Can I get rid of this for loop using apply?

Hey all,

The code below creates a partial dependence plot for the variable x1 in the
linear model y ~ x1 + x1^2 + x2.

I have noticed that the for loop in the code takes a long time to run if the
size of the data is increased.  Is there a way to change the for loop into
an apply statement?  The tricky part is that I need to change the values of
x1 in each step of the loop to give me the appropriate dataset to make
predictions on < cbind(m[,-match("x1",names(m))],x1=a[1,i+1]) >.
If I try
and add the 1112 columns to the dataset a priori and use apply, the code
won't work because the predict function needs the column labeled
"x1".  I
realize I could just grab the form of the linear function and use that
instead of predict(), but I don't want to do that because I want to make
this code applicable to generic model fits.

#create fake data and fit a simple linear regression model
x1 <- rep(c(1,3,4),100)
x2 <- rep(c(1:6),50)
y <- 40+2*x1^.5 - 6*x2 + rnorm(100,0,2)
m <- as.data.frame(cbind(y,x1,x2))
lm1 <- lm(y~x1+I(x1^2)+x2,data=m)

#super small version of R code for partial dependence plot
a <- rbind(c(0:1111)*(max(m$x1)-min(m$x1))/1111 +
min(m$x1),c(0:1111)*0-99999)
for(i in c(0:1111))
{
  a[2,i+1] <-
mean(predict(lm1,cbind(m[,-match("x1",names(m))],x1=a[1,i+1])))
}
plot(a[1,],a[2,],xlab="x1",ylab="Response",type="l",main="Partial
Dependence
Plot")

Many thanks,

Mike

	[[alternative HTML version deleted]]

The answer to my post is yes (which I just figured out).

Solution:



#super small version of R code for pd plot using apply

a <- rbind(c(0:1111)*(max(m$x1)-min(m$x1))/1111 +
min(m$x1),c(0:1111)*0-99999)

b <- matrix(rep(c(0:1111)*(max(m$x1)-min(m$x1))/1111 + min(m$x1), nrow(m)),
nrow(m), 1112, byrow=T)

a[2,] <- apply(b,2,FUN=function(x)
{mean(predict(lm1,cbind(m[,-match("x1",names(m))],x1=x)))  })

plot(a[1,],a[2,],xlab="x1",ylab="Response",type="l",main="Partial
Dependence
Plot")


Mike Dugas

	[[alternative HTML version deleted]]

Great suggestions.  I tested the code on an example and the run time was
reduced from 1 min 12 sec to 3 sec.  Also, I like the suggestion to look at
the quantiles.  I will see what insight it provides in terms of detecting
masked interactions.

I have a couple questions about your code.
First, why not use
xs <- seq(min(x1), max(x1), length = 100) instead of
xs <- with(m, seq(min(x1), max(x1), length = 100)) ?
Second, what is the function geom_line()?  I couldn't find it.

Thanks,
Mike


On 4/23/08, hadley wickham <h.wickham@gmail.com>
wrote:
>
> On Wed, Apr 23, 2008 at 8:33 PM, hadley wickham <h.wickham@gmail.com>
> wrote:
> > > Sure, I am creating a partial dependence plot (reference
Friedman's
> >  > stochastic gradient paper from, I want to say, 2001).  The idea
is to
> find
> >  > the relationship between one of the predictors, say x1, and y by
> creating
> >  > the following plot: take a random sample of actual data points,
hold
> other
> >  > predictors fixed (x2-xp), vary x1 across its range, create a
string of
> >
> >  Put your code doesn't have a random component - you're trying
to
> >  calculate everything combination of the new x_n and the existing
data?
> >   Is that right?
>
> And why are you using so many different values of the x variable?
> 100's should be sufficient to get a smooth curve, not thousands. 
I'd
> also think about displaying not just the mean, but a selection of
> quantiles as well:
>
> Here's one approach:
>
> model <- lm(y ~ poly(x1, 2) + x2, data = m)
>
>
> xs <- with(m, seq(min(x1), max(x1), length = 100))
>
> library(reshape)
> newdf <- expand.grid.df(data.frame(x1 = xs), m[, c("x2"),
drop=F])
>
> predictions <- predict(model, newdata = newdf)
> avg_pred <- tapply(predictions, newdf$x1, mean)
> low_pred <- tapply(predictions, newdf$x1, quantile, 0.25)
> high_pred <- tapply(predictions, newdf$x1, quantile, 0.75)
>
> library(ggplot)
> qplot(xs, avg_pred, min = low_pred, max = high_pred,
geom="ribbon") +
> geom_line()
>
>
> But following your code, it's exhaustive, not random.  This should be
> a little faster because all the predictions are done in one go.
>
> Hadley
>
> --
> http://had.co.nz/
>
	[[alternative HTML version deleted]]