# Failure to unscale coefficients from a linear model with interactions

I have data that require scaling fixed effects to fit linear mixed models with the lmer function from the lme4 package that converge. I want to unscale the estimated coefficients to present results in useful units. Previous posts address this issue, but I am unable to successfully unscale coefficients from models including categorical predictors and interaction terms. I used lmer() with my data set, but I have created an example using lm and the diamonds data for simplicity.

Below is some code that runs without error, but demonstrates that I can't get unscaling to work correctly:

library(ggplot2)
library(dplyr)

# function to rescale variables around 0
zscore <- function(x) (x - mean(x)) / sd(x)

# creates new dataframe with scaled values included
diamonds_s <- diamonds %>% select(carat, cut, price, depth) %>% mutate_each(funs(s = zscore(.)), -cut)


Create unscaled and scaled models

# Linear model with unscaled interaction specified
m1 <- lm(price ~ carat * depth * cut, data= diamonds_s)
b1<- coef(m1)
# summary(m1)

# Linear model with scaled interaction specified
m1s <- lm(price_s ~ carat_s * depth_s * cut, data= diamonds_s)
b1s <- coef(m1s)
# summary(m1s)


Function to rescale coefficients from answer by Ben Bolker

# https://stackoverflow.com/questions/24268031/unscale-and-uncenter-glmer-parameters

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


Use model.matrix to get predictor variables

X <- model.matrix(price ~ carat * depth * cut, data= diamonds_s)

# matrix of response (m) and predictor (m[-1]) means
m <- c(mean(diamonds_s$price), colMeans(X)[-1] ) s <- c(sd(diamonds_s$price), apply(X, 2, sd)[-1] )

# or with unscaled response (not run because I scaled the response, price_s):
# m <- c(0, colMeans(X)[-1] )
# same for standard deviations
# s <- c(1, apply(X, 2, sd)[-1] )


Now attempt to unscale

b1r <- rescale.coefs(b1s, m, s)


The values obtained are wrong, however.

all.equal(b1, b1r, check.names=FALSE, tolerance=0.001)