Standard error of model residual SD from simulation

Question

I'm futzing with getting the SE of model residual standard deviations from a linear regression, and keep getting narrower errors than I should - and I'd like to figure out why.

The basic approach I'm taking is to fit a linear model. Draw simulated coefficients from a multivariate normal distribution. Calculate the RSS and from that, use sqrt(RSS/(n-2)) to calculate the model residual SD. Rinse and repeat 1K times, and then get the SD of the model residual SD.

But... I keep finding that I'm off by an order of magnitude at least. Here's an example in R.

First, the model.

library(palmerpenguins)
plot(bill_length_mm ~ flipper_length_mm, data = penguins)

pen <- lm(bill_length_mm ~ flipper_length_mm, data = penguins)

Then, the simulated coefficients.

library(mvtnorm)
coefTab <- data.frame(rmvnorm(n, coef(pen), sigma = vcov(pen)))

Now, on to getting the residual SD. A function!

get_resid_sd <- function(a, b, y, x){
  pred <- a+b*x
  res <- y - pred
  n <- sum(!is.na(res))
  sqrt(  sum(res^2, na.rm=T)/(n-2) )
}

And let's apply it to our coefficients

coefTab$sigma <- sapply(1:n,
       function(i){
         get_resid_sd(coefTab[i,1], coefTab[i,2], penguins$bill_length_mm,
                      penguins$flipper_length_mm)
       })

Now, the SE

> sd(coefTab$sigma)
[1] 0.01208458

OK, but.... to validate, let's use rstanarm as it naturally produces simulations with no extra work, and can produce equivalent results to lm() using stan_glm() with the appropriate optimizer and null priors.

library(rstanarm)
penStan <- stan_glm(bill_length_mm ~ flipper_length_mm, data = penguins,
                    algorithm = "optimizing", prior = NULL,
                    prior_intercept = NULL, prior_aux = NULL)

Now, the SE of the model residual SD....

> sd(as.matrix(penStan)[,3])
[1] 0.1639105

Huh. You see why I'm a) glad I checked myself and b) why I'm worried I did something tragically wrong.

Would love to know folks' thoughts, as if fixed, I think this is a killer example for students. Feel like I'm falling down on something obvious due to 2020 brain.

Answer 1

$\text{SE}(\hat{\sigma})$ is not the uncertainty in $\hat{\sigma}$ induced by uncertainty in the parameter estimates, which is what you estimated. It's the uncertainty in $\hat{\sigma}$ coming from the data generating process. If we take repeated samples from the population, how does the estimate $\hat{\sigma}$ vary across them?

We can simulate this by repeatedly creating random data, fitting a model to the data, finding the standard deviation of the residuals, and looking at the distribution of those standard deviations.

make_some_data <- function(NOBS = 1000){
  beta <- c(-2,0.4)
  x <- rnorm(n = NOBS, mean = 5, sd = 1)
  y <- beta[1] + beta[2]*x + rnorm(n = NOBS)
  data.frame(x,y)
}
get_model_resids <- function(some_data){
  mod <- lm(y ~ x, data = some_data)
  residuals(mod)
}

estimatedsds <- rep(NA,1000)
for(i in 1:1000){
  estimatedsds[i] <- sd(get_model_resids(make_some_data()))
}
mean(estimatedsds)
# [1] 0.9988542   
sd(estimatedsds)
# [1] 0.02217692

To verify this, we can compare with a stan_glm fit on similar data.

library(rstanarm)
foo <- make_some_data()
stanmod <- stan_glm(y ~ x, data = foo,
                    algorithm = "optimizing", prior = NULL,
                    prior_intercept = NULL, prior_aux = NULL)
sd(as.matrix(stanmod)[,3])
# [1] 0.02255132

We can also compare with a bootstrap estimate of $\text{SE}(\hat{\sigma})$, which is a common way to estimate it.

library(boot)
boot_sds <- function(mydata, myindx){
  mydata <- mydata[myindx,]
  mod <- lm(y ~ x, data = mydata)
  sd(residuals(mod))
}
bootedsds <- boot(foo,boot_sds,R = 10000)
sd(bootedsds$t)
# [1] 0.02239548

You may also be interested in this link, which discusses the exact value of standard error of $\hat{\sigma}$ What is the standard error of the sample standard deviation?