Let $x_i, i=1,\ldots,p$ denote your unscaled predictors. Then the scaled predictors are $z_i = (x_i-m_i)/s_i$, where $m_i$ and $s_i$ are the sample mean and standard deviations of $x_i$, respectively. The linear predictor in the intended model is $$\alpha + \sum_{i=1}^p\beta_i x_i = \alpha + \sum_{i=1}^p\beta_i (s_i z_i + m_i) = (\alpha+\sum_{i=1}^p \beta_i m_i) + \sum_{i=1}^p (s_i\beta_i) z_i$$
Let $\alpha^*$ and $\beta^*_i$ denote the coefficients of the scaled model. From the above we have $$\beta_i^* = s_i\beta_i$$ $$\alpha^* = \alpha+\sum_{i=1}^p \beta_i m_i$$
So if you use the scaled variables as predictors, their coefficients get multiplied by $s_i$ (note that we did not scale $y$, so there is no division by $s_y$), and that operation is easy to "undo". Solving for the original coefficients we get $$\beta_i = \beta_i^*/s_i$$ $$\alpha_i = \alpha_i^* - \sum_{i=1}^p (m_i/s_i)\beta_i^*$$
The original intercept is a linear combination of the coefficients of the model with the scaled predictor. To get its standard error, we have to be careful to take into account the covariance matrix of the coefficient estimates. Using matrix notation, let $C=[1, -m_1/s_1, \ldots,-m_p/s_p]$ be the row vector of weights in the linear combination, and let $V(\hat{b^*})$ be the covariance matrix of the estimates of the model with the scaled predictors $b^*=[\alpha^*, \beta_1^*,\ldots,\beta_p^*]'$. Then $$\hat\alpha = Cb^*, \quad SE(\hat\alpha) = C V(\hat{b^*}) C'$$
In R, you can use, for example the glht function from the multcomp package to do the calculation for you.
Here is an example how this would work with a simple linear model. The only difference with glmer is that you would use fixef instead of coef to extract the model coefficients $\hat{b^*}$, and you would not need to use the df option in glht, because you would use a normal approximation anyway.
> require(mvtnorm)
> set.seed(462627)
> N <- 50
> x <- rmvnorm(N, mean = c(1, -2), sigma = matrix(c(4, 3, 3, 9), nr=2))
> a <- -2; b <- c(0.5, 2)
> y <- rnorm(N, mean= a+x %*% b, sd=1)
>
> # fit model with unscaled predictors
> mod1 <- lm(y ~ x)
> summary(mod1)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-1.44223 -0.60469 0.01115 0.48777 2.73533
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.80351 0.19575 -9.213 4.21e-12 ***
x1 0.35977 0.07559 4.760 1.89e-05 ***
x2 2.03471 0.04348 46.795 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8495 on 47 degrees of freedom
Multiple R-squared: 0.9853, Adjusted R-squared: 0.9846
F-statistic: 1570 on 2 and 47 DF, p-value: < 2.2e-16
>
> # fit model with scaled predictors
> z <- scale(x)
> mod2 <- lm(y ~ z)
> summary(mod2)
Call:
lm(formula = y ~ z)
Residuals:
Min 1Q Median 3Q Max
-1.44223 -0.60469 0.01115 0.48777 2.73533
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.7206 0.1201 -30.97 < 2e-16 ***
z1 0.6573 0.1381 4.76 1.89e-05 ***
z2 6.4619 0.1381 46.80 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8495 on 47 degrees of freedom
Multiple R-squared: 0.9853, Adjusted R-squared: 0.9846
F-statistic: 1570 on 2 and 47 DF, p-value: < 2.2e-16
>
> # calculate intercept from mod1 using mod2
> m <- attr(z, "scaled:center")
> s <- attr(z, "scaled:scale")
> weights <- c(1, -m/s)
> # by hand
> (int <- coef(mod2) %*% weights)
[,1]
[1,] -1.803515
> (se.int <- sqrt(weights %*% vcov(mod2) %*% weights))
[,1]
[1,] 0.1957483
>
> #using glht
> require(multcomp)
> summary(glht(mod2, linfct = rbind(weights), df=mod1$df.residual))
Simultaneous Tests for General Linear Hypotheses
Fit: lm(formula = y ~ z)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
weights == 0 -1.8035 0.1957 -9.213 4.21e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Basic intuition: change of units
Let's consider two models. Let's say $y$ is in units degrees.
Model 1: $x_i$ is denoted in units feet.
$$y_i = a + b x_i + \epsilon_i$$
$b$ is in units degrees per feet.
Model 2: $\tilde{x}_i = 12 x_i$, hence $\tilde{x}_i$ is denoted in units inches.
$$y_i = a + \tilde{b} \tilde{x}_i + \epsilon_i$$
$\tilde{b}$ is in units degrees per inch.
You will find that $\tilde{b} \cdot \frac{12\text{ inches}}{1\text{ foot}} = b$.
You're essentially doing the same thing with your scaling by the standard deviation. It's a change of units and linear models handle changes of units in an intuitively sensible way.
Linear transformation of data leads to linear transformation of coefficient estimates
Linear regression is an linear model. Any linear transformation of the input leads to a clearly defined linear transformation of the estimate.
Let $X$ be the $n$ by $k$ data matrix where we have $n$ observations of $k$ regressors.
Let $\tilde{X} = X A'$ so that $\tilde{\mathbf{x}}_i = A\mathbf{x}$.
Then we have:
\begin{align*}
\tilde{\mathbf{b}} &= (\tilde{X}'\tilde{X})^{-1} \tilde{X}'\mathbf{y} \\
&= \left( A X'X A'\right)^{-1} AX'\mathbf{y}\\
&= A'^{-1} (X'X)^{-1} A^{-1} A X'\mathbf{y}\\
&= A'^{-1} \mathbf{b}
\end{align*}
So if your data is transformed by the linear transformation $A$ so that $\tilde{\mathbf{x}}_i = A \mathbf{x}$ then your estimate $\mathbf{b}$ is transformed so that:
$$\tilde{\mathbf{b}} = A'^{-1} \mathbf{b}$$
Is standardizing (i.e. subtracting mean and scaling by standard deviation) a linear transformation? No if your data does not include a constant, but yes if it does! (i.e. the first column of $X$ is a column of $1$s.)
Example: standardizing with 2 variables and a constant
Let
$$ A = \begin{bmatrix} 1 & 0 & 0 \\ -\frac{\mu_1}{\sigma_1} & \frac{1}{\sigma_1} & 0 \\
-\frac{\mu_2}{\sigma_2} & 0 & \frac{1}{\sigma_2} \end{bmatrix} \quad \quad \mathbf{x} = \begin{bmatrix} 1 \\ x_1 \\ x_2 \end{bmatrix}$$
Then:
$$ \tilde{\mathbf{x}} = A\mathbf{x} = \begin{bmatrix} 1 & 0 & 0 \\ -\frac{\mu_1}{\sigma_1} & \frac{1}{\sigma_1} & 0 \\
-\frac{\mu_2}{\sigma_2} & 0 & \frac{1}{\sigma_2} \end{bmatrix} \begin{bmatrix} 1 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ \frac{x_1 - \mu_1}{\sigma_1} \\ \frac{x_2 - \mu_2}{\sigma_2} \end{bmatrix}$$
That is linear transformation $A$ will standardize your right hand side variables (i.e. subtract mean and divide by standard deviation).
You can show:
$$ A'^{-1} = \begin{bmatrix} 1 & \mu_1 & \mu_2 \\ 0 & \sigma_1 & 0 \\ 0 & 0 &\sigma_2 \end{bmatrix}$$
Hence:
$$\tilde{\mathbf{b}} = A'^{-1} \mathbf{b} = \begin{bmatrix} b_0 + \mu_1 b_1 + \mu_2 b_2 \\ \sigma_1 b_1 \\ \sigma_2 b_2 \end{bmatrix}$$
You get a conceptually similar logic/result if you apply a linear transformation to $y$ (see @whuber's comment for roadmap).
Best Answer
The easiest solution would be to scale your independent variable by dividing it by its standard deviation prior to model fitting. The coefficients and standard errors will then automatically be on the standardized scale.
i) You should definitely report the standard error as a measure of uncertainty. Additionally you can report the standard deviations in a descriptive analysis of the independent variables
ii) Proceed as I said above, or use the delta method, i.e. multiply the standard errors obtained through the "traditional method" by the inverse of the square of the standard deviation.