Conceptually, what is the bias of the standard error of an estimator?

Question

I'm reading Muthén and Muthén (2002) to learn how to use Monte Carlo simulation to estimate statistical power in regards to the coefficients of a model that is linear in its coefficients.

I understand the bias of an estimator is the difference between its expected value versus the population parameter's value, i.e., $$\huge\mathrm{E}[\hat\theta]-\theta$$

However, I don't think I conceptually understand what is the bias of the standard error of an estimator, i.e., the bias of $\sigma_{\hat\theta}$.

Here is the excerpt from page 606 of Muthén and Muthén (2002)$^{\text{see footnote}}$:

Parameter bias is evaluated using the information in columns 1 and 2 of the output. The column labeled starting gives the population parameter values. The column labeled average gives the parameter estimate average over the replications of the Monte Carlo study. For example, the first number in column 2, 0.7963, is the average of the factor loading estimates for y1 over 10,000 replications. To determine its bias, subtract the population value of 0.8 from this number and divide it by the population value of 0.8. This results in a bias of –0.005, which is negligible.

Standard error bias is evaluated using the information in columns 3 and 4 of the output. The column labeled SD gives the standard deviation of each parameter estimate over the replications of the Monte Carlo study. This is considered to be the population standard error when the number of replications is large. The column labeled SE Average gives the average of the estimated standard errors for each parameter estimate over the replications of the Monte Carlo study. Standard error bias is calculated in the same way as parameter estimate bias as described previously.

Thus, is it the difference between the sample-derived $\huge\sigma_{\hat\theta}$ versus what $\huge\sigma_{\hat\theta}$ would have been if the estimator could hypothetically be applied across the entire population? Expressed mathematically with an obnoxiously long subscript:

$$\huge\sigma_{\hat\theta}-\sigma_{\hat\theta_{\text{if }\hat\theta\text{ were to be applied over the entire population}}}$$

Footnote: They divide the bias of the estimator by the parameter value to obtain a useful proportion, but they refer to this quotient as the bias of the estimator, which I believe is inconsistent with how most statisticians define bias, but that's beside the point of my question.

References:

• Muthén, L. K., & Muthén, B. O. (2002). How to Use a Monte Carlo Study to Decide on Sample Size and Determine Power. Structural Equation Modeling: A Multidisciplinary Journal, 9(4), 599–620. https://doi.org/10.1207/S15328007SEM0904_8

Answer 1

The standard error of a parameter estimate is just a model-based estimate of the true SD of the estimator. Often, standard errors are based on asymptotic results that may not be valid for finite samples. The process you describe is a way of checking whether the standard error is unbiased for finite samples.

Specifically, we can generate a huge number of replicates (datasets) and compute the parameter estimate of interest based on each. We can then find the sample SD of these (iid) estimates, which is consistent for the true SD of the estimator (i.e., as you generate more replicates, the sample SD will approach the true SD with probability 1). So, with enough replicates, you can think of this sample SD as the true SD.

Now, with each replicate, you also get a SE of the parameter estimate of interest (an estimate of the true SD of the estimator). The average of these SEs will approach the mean SE as the number of replicates grows (with probability 1). If the estimator of the true SD is unbiased, as the number of replicates grows, the average SE will equal the true SD (which will equal the sample SD described above).

Essentially, this process allows us to determine whether the model-based SE is a good approximation to the true SD of the parameter estimator.

Answer 2

The simulation is like a bootstrap from the population, so the sampling variation across the bootstrap simulations is a way to calculate the standard error of a population of that size.

The standard deviation of the estimatr should be approximately equal to the mean standard error.

You can demonstrate this with something less complex than SEM.

Using R, generate a sample of n = 10 from a normal distribution.

library(dplyr)

n <- 10

res <- lapply(1:10000, function(i) {
  x <- rnorm(n)
  se <- sd(x) / sqrt(n)
  m <- mean(x)
  return(c(m = m, se = se))
}) 

df <- dplyr::bind_rows(res)

mean(df$se)
sd(df$m)

The mean standard error and the standard deviation of the means are very similar:

> mean(df$se)
[1] 0.3072158
> sd(df$m)
[1] 0.3162565

I can conclude that the estimate of the standard error is (about) right.

We didn't even need to use the standard deviation, because I know the population standard deviation that I used (it's 1.0).

> 1/sqrt(10)
[1] 0.3162278

In something like CFA, I don't know the population variance of the parameter.

Generate them from a non-normal distribution, say a $\chi^2_1$:

res <- lapply(1:10000, function(i) {
  x <- rchisq(n, 1)
  se <- sd(x) / sqrt(n)
  m <- mean(x)
  return(c(m = m, se = se))
}) 

df <- dplyr::bind_rows(res)

mean(df$se)
sd(df$m)

Now my sd and mean(se) are different:

> mean(df$se)
[1] 0.3928587
> sd(df$m)
[1] 0.4440972

Do it again with a sample of 100, instead of 10:

> mean(df$se)
[1] 0.1393581
> sd(df$m)
[1] 0.1429515

The values are much more similar to each other. With a small sample and a non-normal distribution, the regular formula for the standard error is biased. When the sample size increases, the bias is reduced.