How do regression results change after standardization, as a general rule? Based on the simulation below, it appears that standardizing all variables in a data set affects OLS results in the following ways:   


*

*Coefficient estimates change

*Standard errors change

*P-values remain the same except the p-value for the intercept coefficient.


Are these results general? Do they also apply to other models, not only OLS? 
Simulation:
#Standardization simulation
remove(list = ls())

set.seed(42)
n = 50
t <- rnorm(n, mean = 2.2, sd = 5)
x1 <- rnorm(n, mean = 1.5, sd = 5)
x2 <- rnorm(n, mean = 3.3, sd = 6)
x3 <- rnorm(n, mean = 2, sd = 7)

betas <- matrix(runif(4, min = -5, max = 5))

inputs <- as.matrix(cbind(t, x1, x2, x3))

y <- (inputs %*% betas) + rnorm(n, mean = 0, sd = 20)

data <- data.frame(cbind(y, inputs))

standardize <- function(variable){
    demeaned <- variable - mean(na.omit(variable))
    sd <- sqrt(var(na.omit(variable)))
    return(demeaned/sd)
}

stan.data <- data.frame(apply(data, 2, FUN = standardize))

summary(lm(y ~ t + x1 + x2 + x3, data = data))

summary(lm(y ~ t + x1 + x2 + x3, data = stan.data))

 A: This method is called beta coefficients (at least by Wooldridge). We compute the Z-scores of all variables, and then run the regression on the transformed data. Now what happens?
The slopes are different: before we had a slope of say $\hat \beta_1$, but now we have a slope of $\hat \sigma_1 /\hat \sigma_y \cdot \hat \beta_1$, where $\hat \sigma_i$ is the sample standard deviation. 
There is no intercept; or rather it will be estimated as very (very) close to 0. Hence the change in p-value of the intercept you mentioned.
The interpretation of $\hat \beta_i$ is now: what happens with y, when we increase the standard deviation of $x_i$ by 1. In a sense, we can say something about the most "important" variable (but I dislike that word).
One final note: just as you you stated, the p-values are not going to be different and so you will always have the same significance as before. 
A: Yes, they are general for OLS.
Let me take the example of a simple regression. We know the slope coefficient is
$$
\hat\beta_1=\frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})}{\sum_i(x_i-\bar{x})^2}
$$
or
$$
\hat\beta_1=\frac{\widehat{Cov}(y_i,x_i)}{\widehat{Var}(x_i)}
$$
and the intercept estimate is
$$
\hat\beta_0=\bar y-\hat\beta_1\bar x.
$$
Standardizing yields $\bar x=\bar y=0$ and $\widehat{Var}(x_i)=\widehat{Var}(y_i)=1$. Hence, 
$$
\hat\beta_1=\widehat{Cov}(y_i,x_i)=\widehat{Corr}(y_i,x_i)
$$
and 
$$
\hat\beta_0=0.
$$
The standard errors are the square roots of the diagonal elements of $s^2(X'X)^{-1}$, e.g., $s.e.(\hat\beta_1)=s/(\sum_i(x_i-\bar{x})^2)^{1/2}$.
In the scaled case, this standard error then becomes
$$\tilde s/(n-1)^{1/2},$$
where $\tilde s$ is the square root of the error variance estimate of the scaled regression. We get $n-1$ in the denominator because $\widehat{Var}(x_i)=1/(n-1)\sum_i(x_i-\bar{x})^2=1$, so $\sum_i(x_i-\bar{x})^2=n-1$. Hence, the standard errors are different.
The $t$-ratios are nevertheless the same. For the unscaled case, write
\begin{align*}
t&=\frac{\hat{\beta_1}}{s.e.(\hat\beta_1)}\\
&=\frac{\frac{\widehat{Cov}(y_i,x_i)}{\widehat{Var}(x_i)}}{s/(\sum_i(x_i-\bar{x})^2)^{1/2}}\\
&=\frac{\frac{\widehat{Cov}(y_i,x_i)}{\widehat{Var}(x_i)}}{s/((n-1)^{1/2}sd(x))}\\
&=(n-1)^{1/2}\frac{\widehat{Cov}(y_i,x_i)sd(x)}{s\widehat{Var}(x_i)}\\
&=(n-1)^{1/2}\frac{\widehat{Corr}(y_i,x_i)\widehat{Var}(x_i)sd(y)}{s\widehat{Var}(x_i)}\\
&=(n-1)^{1/2}\frac{\widehat{Corr}(y_i,x_i)sd(y)}{s},
\end{align*}
where the 5th equality follows from rearranging 
$$
\widehat{Corr}(y_i,x_i)=\frac{\widehat{Cov}(x_i,y_i)}{sd(x)sd(y)}
$$
For the unscaled case, write
\begin{align*}
\tilde t&=\frac{\widehat{Corr}(y_i,x_i)}{\frac{\tilde s}{(n-1)^{1/2}}}\\
\end{align*}
An application of the Frisch Waugh Lovell theorem will demonstrate that scaling the regressor does nothing to the residuals. Scaling $y$, however, gives residuals $\tilde u_i$ which are related to the unscaled ones $\hat u_i$ via $\tilde u_i=\hat u_i/sd(y)$. Hence, $\tilde s=s/sd(y)$. Thus,
\begin{align*}
\tilde t&=\frac{(n-1)^{1/2}\widehat{Corr}(y_i,x_i)}{\frac{s}{sd(y)}}\\
&=\frac{(n-1)^{1/2}\widehat{Corr}(y_i,x_i)sd(y)}{s}=t\\
\end{align*}
If the test statistics are the same, so will of course the p-values.
That the standard error on the intercept is different follows from the fact that the scaled intercept is estimated to be zero very precisely.
P.S.: scale achieves the goals of your function, too.
