Failure to unscale coefficients from a linear model with interactions

Question

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.

related threads: https://stackoverflow.com/questions/23642111/how-to-unscale-the-coefficients-from-an-lmer-model-fitted-with-a-scaled-respon

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

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[1]*beta[-1]/sigma[-1]
    beta2[1] <- sigma[1]*beta[1]+mu[1]-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[1]) 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)

I appreciate any input about this issue.

Answer 1

I think I had the same problem. I've solved scaling the interactions terms using their own means and standard deviations.

This thread may be useful, check the detailed Ben Bolker answer.