How to interpret coefficients that flip sign once standardized?

Question

Take the simple simulation below.

set.seed(3424131)
x_1 = rnorm(100, mean = 10, sd = 4)
x_2 = rnorm(100, mean = 3, sd = 2)


covariate = 0.05-0.03*x+0.04*x_2

y = rnbinom(n= length(x_1), mu = exp(covariate), size = 50)

temp = glm.nb(y~x_1+x_2+x_1*x_2)

summary(temp)


x_new_1 = scale(x_1)
x_new_2 = scale(x_2)

temp_2 = glm.nb(y~x_new_1+x_new_2+x_new_1*x_new_2)

summary(temp_2)

Before standardizing(temp_1), the coefficients are as follow:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.01450    0.60287  -1.683   0.0924 .
x_1          0.05332    0.05079   1.050   0.2937  
x_2          0.26133    0.13552   1.928   0.0538 .
x_1:x_2     -0.01718    0.01235  -1.391   0.1643  

And after standardizing(temp_2), the coefficients are as follow.

Coefficients:
                Estimate Std. Error z value Pr(>|z|)
(Intercept)     -0.18402    0.11300  -1.629    0.103
x_new_1         -0.01949    0.11293  -0.173    0.863
x_new_2          0.14715    0.10979   1.340    0.180
x_new_1:x_new_2 -0.14799    0.10640  -1.391    0.164

After standardizing, the sign of the coefficient for $x_1$ flipped. How do we interpret this?

It seems like before standardizing, an increase in $x_1$ resulted in the increase of the mean, but after standardizing, an increase in $x_1$ results in a decrease of the mean, which seems paradoxical.

Also, does the change of sign have to do with the interaction term? I couldn't replicate the flipping of the sign without the interaction term.

Answer 1

Your first model supposes the parameters of the response distribution are a linear combination of four variables: an intercept, two variables $x_i,$ and their "interaction" $x_1x_2.$ Such a linear combination takes the form

$$\beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_{12}x_1x_2.$$

Your second model re-expresses the $x_i$ in terms of their standardized forms, $x_i = \sigma_i z_i + \mu_i,$ which when substituted into the first model give

$$\beta_0 + \beta_1 (\sigma_1 z_1 + \mu_1) + \beta_2 (\sigma_2 z_2 + \mu_2) + \beta_{12} (\sigma_1 z_1 + \mu_1) (\sigma_2 z_2 + \mu_2) .$$

When you use the usual rules of algebra with real numbers to express this as a linear combination of the new variables $\alpha_0+\alpha_1z_1+\alpha_2z_2+\alpha_{12}z_1z_2$ you find

$$\begin{aligned} \alpha_0 &= \beta_0 + \beta_1\mu_1 + \beta_2 \mu_2 + \beta_{12}\mu_1\mu_2;\\ \alpha_1 &= \beta_1\sigma_1+\beta_{12}\sigma_1\mu_2;\\ \alpha_2 &= \beta_2\sigma_2 + \beta_{12}\sigma_2\mu_1;\text{ and}\\ \alpha_{12} &= \beta_{12}\sigma_1\sigma_2. \end{aligned}$$

You have noticed that the signs of $\beta_1$ and $\alpha_1$ are opposite. Although multiplying $\beta_1$ by $\sigma_1$ won't change its sign ($\sigma_1$ is the standard deviation of the $x_1$ values), adding a multiple of $\beta_{12}\mu_2$ definitely can do this, depending on the signs and magnitudes of (a) the original estimated interaction term and (b) the mean of $z_2.$