Interpreting the coefficients when predictors are standarised Let's say I have a model of the form:
mdl <- lm(y ~ x1 + x2 + x3, data=dat)

I am interested in knowing the relationship between x3 and y when x1 and x2 remain unchanged. To do this, I do the following:
dat$x1 <- mean(dat$x1)
dat$x2 <- mean(dat$x2)
dat$x3 <- dat$x3

pred <- predict(mdl, data=dat)
mean(pred) 

If I interpret mean(pred) correctly, it represents average change in y with a unit change in x3 while keeping x2 and x1 constant. Please correct me if I am wrong.
Now imagine a second model wherex1, x2 and x3 are standardised.
dat$z.x1 <- scale(dat$x1,scale = T,center = T)
dat$z.x2 <- scale(dat$x2,scale = T,center = T)
dat$z.x3 <- scale(dat$x3,scale = T,center = T)

mdl1 <- lm(y ~ z.x1 + z.x2 + z.x3, data=dat)

# keeping x1 and x2 at their mean values (which is zero since the variables 
are standardised)  

dat$z.x1 <- 0
dat$z.x2 <- 0
dat$z.x3 <- dat$z.x3

pred.z <- predict(mdl1, data=dat)

mean(pred.z) 

This is where my question is. Does mean(pred.z) has the same interpretation as mean(pred) in the previous model? Can mean(pred.z) be interpreted as mean change in y with a unit change in x3 while keeping x2 and x1 constant. OR 
Does mean(pred.z) represent mean change in y when x3 changes by one standard deviation since I had standardised x3 in this model? 
Thank you
 A: You are doing more work than necessary to know "the relationship between x3 and y when x1 and x2 remain unchanged," and at first glance you seem to be making an error in your approach (if I understand it properly).
The regression coefficient for x3 in your first model is precisely the change in y per unit change in x3 with the other predictors held constant. In your second model, with x3 standardized before the regression, the regression coefficient for z.x3 is the change in y per one-standard-deviation change in x3, again with the other predictors held constant.
Using predict() and mean():
pred <- predict(mdl, data=dat)
mean(pred) 

after setting all the other predictors to their mean values doesn't seem, at first glance, to give you what you are looking for. It gives you the mean y prediction for a data set that kept the original x3 values but kept all of the other predictors at their individual means. The distribution of predicted y values would thus be very dependent on the distribution x3 values and wouldn't necessarily bear any relation to "the relationship between x3 and y when x1 and x2 remain unchanged." It's easy to come up with a scenario in which the mean of predicted y values would be 0 even with a strong linear relation between y and x3. If I'm missing something in your approach, please clarify.
In your second model, the pre-scaling of x3 will lead to a corresponding change in the regression coefficient for z.x3 so it seems that the final predicted y values with 
pred.z <- predict(mdl1, data=dat)
mean(pred.z)

should be the same numerically as in the first model, although I haven't checked that directly.
Again, the regression coefficients themselves provide directly your answer.
