transformaciones de coeficientes en regresiones

Cómo hacemos para volver a la escala original una vez que hicimos una regresión con variables centradas y estandarizadas?

ejemplo:

set.seed(36)
n = 1000
a_origi = runif(1, -2,2)
b_origi = runif(1, -2,2)
sd = runif(1)
x = runif(n, 0, 5)
y = rnorm(n, a_origi + b_origi * x, sd )
plot(x, y)
abline(a_origi, b_origi)

Hacemos una regresión con la variable \(x\) centrada y estandarizada (como debe ser):

xs = scale(x)
fit = lm(y ~ xs)
print(fit)

## 
## Call:
## lm(formula = y ~ xs)
## 
## Coefficients:
## (Intercept)           xs  
##      2.2359       0.9795

Cómo hacemos para pasar las estimaciones de los coeficientes a la escala original?

Si lo único que nos importa es graficar los resultados, podemos usar este truco:

plot(x, y)
lines(x, fit$coefficients[1] + fit$coefficients[2] * xs)

La pendiente está ahora en unidades de desvío estándar, por lo que para recuperar la escala original tenemos que hacer

b_est = fit$coefficients[2] / sd(x)

La ordenada al origen está en la media de \(x\)

a_est = fit$coefficients[1] - mean(x) * b_est

graficamos

plot(x, y)
abline(a_est, b_est, col = 2)

Hagamos esto muchas veces para ver si funciona

n = 1000
nreps = 5000
adif = numeric(nreps)
bdif = numeric(nreps)
for(i in 1:nreps){
  a_origi = runif(1, -2,2)
  b_origi = runif(1, -2,2)
  sd = runif(1)
  x = runif(n, 0, 5)
  y = rnorm(n, a_origi + b_origi * x, sd )
  xs = scale(x)
  fit = lm(y ~ xs)
  b_est = fit$coefficients[2] / sd(x)
  a_est = fit$coefficients[1] - mean(x) * b_est
  adif[i] = a_est - a_origi
  bdif[i] = b_est - b_origi
}

hist(adif, 100)

hist(bdif, 100)