Standard Error – Is ‘Std. Error’ Provided by lmer in lme4 Actually Standard Deviation?

lme4-nlmemixed modelrstandard deviationstandard error

I have been working with lmer in lme4, and wanted to construct a density curve of the slopes produced by my model (rather than just confidence intervals). The model output of lmer produces a slope, and something referred to as 'std. error'. If I calculate the standard deviation from the standard error, the standard deviation is much higher than I would have expected.

I investigated this further by using confint to select 1000 random confidence intervals, and comparing the outputted slopes to an rnorm using both what I calculated the standard deviation would be, and the standard error given by lmer. As you can see, the standard error is much closer to the slopes from random confidence intervals than it is to the standard deviation. Am I right in thinking that what is referred to as standard error in the lmer output is, at least for the purposes of the density distribution of slopes, standard deviation?

Some example code:

library(lme4)

set.seed(123)

#generate random dataset
df <- data.frame(weight=(rnorm(10000, 1, 0.3)),
                 height=(rnorm(10000, 5, 6)),
                 species = rep(c("mudkip", "torchic", "treecko"), length.out = 10000)) 
# model

mod <- lmer(data=df, weight ~ height + (1|species), REML=F)

slope <- coef(summary(mod))[ , "Estimate"][2]
se <- coef(summary(mod))[ , "Std. Error"][2]

sd <- se * sqrt(10000)

# Randomly extract slopes for varying confidence intervals

slope_num <- 1
ci_slope <- c()
for(i in runif(1000)){
  ci <- confint(mod,method="Wald", level=i)
  lwr <- colnames(ci)[1]
  upr <- colnames(ci)[2]
  ci_slope_lwr <- ci["height", lwr]
  ci_slope_upr <- ci["height", upr]
  slope_num <- slope_num + 1
  ci_slope <- c(ci_slope, ci_slope_lwr, ci_slope_upr)
}

# Generate rnorms using Std. Error and calculated standard deviation

rnorm_se <- rnorm(1000, slope, se)
rnorm_sd <- rnorm(1000, slope, sd)

# Comparative plots

plot(density(rnorm_se), lwd=2)
lines(density(ci_slope), col="blue")
lines(density(rnorm_sd), col="red") # I can't construct a readable graph that can fit sd with anything else

Plot of rnorm using 'Std. Error', rnorm using standard deviation, and density of slopes from random confidence intervals

Best Answer

The "standard error" of a parameter is the standard deviation of the parameter's sampling distribution (e.g. wiki). In other words, theoretically, if you were to sample new data and rerun the analysis again and again, you would end up with a distribution of slope parameters, and standard error is an estimate of the standard deviation of that distribution.

Where I think you're confused is that when the parameter in question is the sample mean, there is a simple relationship between the sample standard deviation and the standard error of the mean (that is, the standard deviation of the sampling distribution of the estimated mean): $\text{SEM} = \frac{SD}{\sqrt{n}}$. This relationship doesn't apply to parameter estimates from regression models.

Related Solutions

Regression – How to Calculate the Specific Standard Error Relevant for a Specific Point Estimate Within a Linear Regression

The first part of my question, is how do you calculate this specific Standard Error at a specific point estimate?

You don't specify if you mean simple linear or multiple regression. I'll assume the general case. Let's do it at a point $x^* = (1,x_1^*,x_2^*,...,x_p^*)$

$$\text{Var}(\hat y^*) = \text{Var}(x^*\hat\beta)= \text{Var}(x^*(X^TX)^{-1}X^T y)$$

$$= x^*(X^TX)^{-1}X^T \text{Var}(y) X(X^TX)^{-1}x^{*T}$$ $$ = \sigma^2 x^*(X^TX)^{-1}X' I X(X^TX)^{-1}x^{*T} $$ $$= \sigma^2 x^*(X^TX)^{-1}x^{*T}$$

If $h^*_{ii} = [x^*(X^TX)^{-1}x^{*T}]_{ii}$ then $\text{Var}(\hat y_i) = \sigma^2 h^*_{ii}$.

Of course, $\sigma^2$ is unknown and must be estimated.

The standard error is the square root of that estimated variance up above.

Could one provide a link to a numerical example to facilitate my interpretation of the formula?

I'll try to dig one up.

My second part to this overall question is: How come the resulting hourglass shape of the resulting Confidence Interval as depicted does not break the linear regression assumption that the variance of residuals remain constant across observations (the heteroskedasticity thing)?

1) it's a confidence interval for where the mean is, not the variance of the data; it reflects our uncertainty in the parameters as they feed through (via the design, $X$) to the the estimate of the mean. Something assumed true for one thing not being true for a different thing doesn't violate the assumption for the first thing.

2) Your statement "the linear regression assumption that the variance of residuals remain constant across observations" is factually incorrect (though I know what you're getting at). That is not an assumption of regression - in fact, outside specific cases, it's untrue for regression. What is assumed constant is the variance of the unobserved errors. The variance of the residuals is not constant. In fact it 'bows in' in opposite fashion to the way the variance above 'bows out', both due to the phenomenon of leverage.

Edits in response to followup questions:

Why would the variance bow in?

I'll do it algebraically and then expand on the explanation in the text above:

\begin{eqnarray} \text{Var}(y-\hat y) &=& \text{Var}(y) + \text{Var}(\hat y) - 2 \text{Cov}(y,\hat y)\\ &=&\sigma^2 I + \text{Var}(X \hat \beta) - 2 \text{Cov}(y,X \hat \beta)\\ &=&\sigma^2 I + \text{Var}(X (X^TX)^{-1}X^T y) - 2 \text{Cov}(y,X (X^TX)^{-1}X^T y)\\ &=&\sigma^2 I + X (X^TX)^{-1}X^T\text{Var}(y) X (X^TX)^{-1}X^T - 2 \text{Cov}(y, y)X (X^TX)^{-1}X^T\\ &=&\sigma^2 I + X (X^TX)^{-1}X^T(\sigma^2 I) X (X^TX)^{-1}X^T - 2 \sigma^2 I X (X^TX)^{-1}X^T\\ &=&\sigma^2 I + \sigma^2 X (X^TX)^{-1}X^T X (X^TX)^{-1}X^T - 2 \sigma^2 I X (X^TX)^{-1}X^T\\ &=&\sigma^2 I + \sigma^2 X (X^TX)^{-1}X^T - 2 \sigma^2 X (X^TX)^{-1}X^T\\ &=&\sigma^2 [I + X (X^TX)^{-1}X^T - 2 X (X^TX)^{-1}X^T]\\ &=&\sigma^2 [I - X (X^TX)^{-1}X^T]\\ &=& \sigma^2(I-H) \end{eqnarray}

where $H = X(X^TX)^{-1}X^T$. Therefore the variance of the $i^{\tt{th}}$ residual is $\sigma^2(1-h_{ii})$ where $h_{ii}$ is $H(i,i)$ (some texts will write that as $h_i$ instead).

As you see, it's smaller, when $h$ is larger... which is when the pattern of $x$-values is further from the center of the $x$'s. In simple regression $h$ is larger when $(x-\bar x)$ is larger.

Now as to why, note that $\hat y = Hy$ ($H$ is called the hat-matrix for this reason).

That is, the fit at the $i^{\tt{th}}$ observation responds to movements in $y_i$ in proportion to $h_{ii}$, or $\frac{\partial \hat{y}_i}{\partial y_i} = h_{ii}$. So when $h$ is larger, $y$ pulls the line more toward itself, making its residual smaller.

There's a more intuitive discussion in the context of simple linear regression here that may help motivate it for you.

I interpret that as large errors near the Mean with smaller errors away from the Mean.

No, we're not discussing errors, they have constant variance. We're discussing residuals. They are not the errors and don't have constant variance; they're related but different.

The bit of material I have read on the subject, suggests just the opposite...

Can you point me to something that does this? Recall that we're discussing the residual variability here.

Additionally, how would you define heteroskedasticity?

Having non-constant variance. That is, when the regression assumption about the variance being constant doesn't hold, you have heteroskedasticity.

See Wikipedia: http://en.wikipedia.org/wiki/Heteroscedasticity

And, what do you mean by the variance of unobserved errors?

You don't observe the errors, but the model assumes they have constant variance, $\sigma^2$. The "variance of unobserved errors" is thus simply "$\sigma^2$".

How can you measure those since they are unobserved?

Individually, you can't, at least not very well. You can roughly estimate them by the residuals, but they don't even have the same variance, as we saw. However, you can estimate their variance reasonably well from the residuals, if you appropriately adjust for the fact that the residuals are on average smaller than the errors.

Solved – How to interpret slope parameter estimates for linear models in R

I think it's always difficult to interpret the size of interaction effects, but it can be done in a different way to make it easier.

If you create two new variables where the first is 0 for diet B and equal to time for diet A, and the other variable is the oppoosite, you assess the effect of the slope for each diet more directly:

TimeA <- ifelse(Diet=="A", Time, 0)
TimeB <- ifelse(Diet=="B", Time, 0)

Nitrogen ~ Diet + TimeA + TimeB + (1|replicate)

You will now get separate estimates for the slope for time A and time B, with standard errors so that you can calculate confidence intervals to report. You will see that the estimates for Diet are the same, and the estimate for TimeA is the same as the estimate for Time in your model. The difference is that in this model, you will get one estimate for the slope of time for Diet B, and it will probably be non-significant in your model.

Note that this assumes that the baseline for Time is 0. If not, the results will be more difficult to interpret.

Best Answer

Related Solutions

Regression – How to Calculate the Specific Standard Error Relevant for a Specific Point Estimate Within a Linear Regression

Solved – How to interpret slope parameter estimates for linear models in R

Related Question