How to unscale the coefficients from an lmer()-mod

I fitted a model in R with the lmer()-function from the lme4 package. I scaled the dependent variable:

    mod <- lmer(scale(Y)
                ~ X
                + (X | Z),
                data = df,
                REML = FALSE)

I look at the fixed-effect coefficients with fixef(mod):

    > fixef(mod)
    (Intercept)      X1          X2         X3           X4 
     0.08577525 -0.16450047 -0.15040043 -0.25380073  0.02350007

It is quite easy to calculate the means by hand from the fixed-effects coefficients. However, I want them to be unscaled and I am unsure how to do this exactly. I am aware that scaling means substracting the mean from every Y and deviding by the standard deviation. But both, mean and standard deviation, were calculated from the original data. Can I simply reverse this process after I fitted an lmer()-model by using the mean and standard deviation of the original data?

Thanks for any help!

Update: The way I presented the model above seems to imply that the dependent variable is scaled by taking the mean over all responses and dividing by the standard deviation of all the responses. Usually, it is done differently. Rather than taking the overall mean and standard deviation the responses are standardized per subject by using the mean and standard deviation of the responses of that subject. (This is odd in an lmer() I think as the random intercept should take care of that... Not to mention the fact that we are talking about calculating means on an ordinal scale...) The problem however stays the same: Once I fitted such a model, is there a clean way to rescale the coefficients of the fitted model?

Updated: generalized to allow for scaling of the response as well as the predictors.

Here's a fairly crude implementation.

If our original (unscaled) regression is

Y = b0 + b1*x1 + b2*x2 ...

Then our scaled regression is

(Y0-mu0)/s0 = b0' + (b1'*(1/s1*(x1-mu1))) + b2'*(1/s2*(x2-mu2))+ ...

This is equivalent to

Y0 = mu0 + s0((b0'-b1'/s1*mu1-b2'/s2*mu2 + ...) + b1'/s1*x1 + b2'/s2*x2 + ...)

So bi = s0*bi'/si for i>0 and

b0 = s0*b0'+mu0-sum(bi*mui)

Implement this:

 rescale.coefs <- function(beta,mu,sigma) {
    beta2 <- beta ## inherit names etc.
    beta2[-1] <- sigma[1]*beta[-1]/sigma[-1]
    beta2[1]  <- sigma[1]*beta[1]+mu[1]-sum(beta2[-1]*mu[-1])
    beta2
 }

Try it out for a linear model:

m1 <- lm(Illiteracy~.,as.data.frame(state.x77))
b1 <- coef(m1)

Make a scaled version of the data:

ss <- scale(state.x77)

Scaled coefficients:

m1S <- update(m1,data=as.data.frame(ss))
b1S <- coef(m1S)

Now try out rescaling:

icol <- which(colnames(state.x77)=="Illiteracy")
p.order <- c(icol,(1:ncol(state.x77))[-icol])
m <- colMeans(state.x77)[p.order]
s <- apply(state.x77,2,sd)[p.order]
all.equal(b1,rescale.coefs(b1S,m,s))  ## TRUE

This assumes that both the response and the predictors are scaled.

If you scale only the response and not the predictors, then you should submit (c(mean(response),rep(0,...)) for m and c(sd(response),rep(1,...)) for s (i.e., m and s are the values by which the variables were shifted and scaled).
If you scale only the predictors and not the response, then submit c(0,mean(predictors)) for m and c(1,sd(predictors)) for s.