The other day in class, while talking about instances (e.g., analyzing clustered data or heteroskedastic residuals) where adjustments are required to the standard errors of a regression model, a student asked: how do we know what the ‘true’ standard error should be in the first place– which is necessary to know if it is too high or too low.

This short simulation illustrates that, over repeated sampling from a specified population, the **standard deviaton** of the regression coefficients can be used as the true standard errors. Take for example the case where the residuals of a regression are heteroskedastic instead of the assumed homoskedastic residuals. Parameter estimates will be unbiased but the standard errors will be biased.

To illustrate this, I create one dataset which shows heteroskedastic residuals:

```
set.seed(2468)
ns <- 400 #sample size
B1 <- 3 #regression coefficient
x1 <- runif(ns, 10, 30)
err <- rnorm(ns, 0, x1^2) #variance [or sd] increases as x1 increases
y <- 50 + B1 * x1 + err
m1 <- lm(y ~ x1) #run a basic regression
summary(m1) #we see that the the coefficient for x1
```

```
##
## Call:
## lm(formula = y ~ x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1980.99 -274.52 7.57 270.21 1458.05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.230 95.180 0.423 0.673
## x1 2.907 4.667 0.623 0.534
##
## Residual standard error: 506.3 on 398 degrees of freedom
## Multiple R-squared: 0.000974, Adjusted R-squared: -0.001536
## F-statistic: 0.388 on 1 and 398 DF, p-value: 0.5337
```

```
plot(x1, rstudent(m1), ylab = 'Studentized residuals') #classic fan spread of residuals
library(lmtest)
```

`bptest(m1) #Breusch Pagan test: results show statistical significance indicating a violation`

```
##
## studentized Breusch-Pagan test
##
## data: m1
## BP = 98.226, df = 1, p-value < 2.2e-16
```

However, the question then becomes: what should the standard error be? So, we can run a simulation where we run a regression 1,000 times (1,000 replications) and collect the output. Typically, using heteroskedasticity corrected standard errors should alleviate the problem which we will see.

```
set.seed(123)
ns <- 400 #sample size
B1 <- 3 #regression coefficient
reps <- 500 #replications
#container to save the coefficient and the standard errors
cfs <- ses <- ses2 <- matrix(NA, ncol = 2, nrow = reps)
library(sandwich) #to get the adjusted standard errors
set.seed(3210) #for reproducability
#the actual simulation
for (i in 1:reps){
x1 <- runif(ns, 10, 30)
err <- rnorm(ns, 0, x1^2) #variance [or sd] increases as x1 increases
y <- 50 + B1 * x1 + err
m2 <- lm(y ~ x1) #run a basic regression
cfs[i, ] <- coef(m2)
ses[i, ] <- sqrt(diag(vcov(m2)))
ses2[i, ] <- sqrt(diag(vcovHC(m2) ))
#cat(i, "::")
}
```

After running the simulation, we can compare the actual standard errors to the empirically determined standard errors. The point estimate is ~ 3.00 (i.e., \(B = 3.04\)).

```
### these are for the point estimates: should be unbiased
colMeans(cfs) #should be around 50 and 3.00
```

`## [1] 47.974758 3.044893`

```
hist(cfs[,2], main = 'Parameter estimates [3.00]')
abline(v = 3, col = 'red')
```

For the standard errors, we can compare the 1) empirical standard errors; 2) the unadjusted standard errors; and 3) the corrected standard errors.

```
### these are for the standard errors
apply(cfs, 2, sd) #these show the empirical standard errors for the intercept and x1
```

`## [1] 78.909899 4.819149`

`colMeans(ses) #these are the model based standard errors`

`## [1] 88.474339 4.247855`

`colMeans(ses2) #these are adjusted standard errors`

`## [1] 78.28279 4.78608`

We can see that the *unadjusted* standard errors for \(x_{1}\) are **too low** (`(4.24 - 4.82)/4.82 = -12`

%) while the *adjusted* standard errors are close to the empirical standard errors (i.e., 4.79 vs. 4.82).

```
dat <- data.frame(se = c(ses[,2], ses2[,2]), gr = rep(c('unadj', 'hc'), each = reps))
library(ggplot2)
ggplot(dat, aes(x = se, group = gr, col = factor(gr))) +
geom_density() +
geom_vline(xintercept = sd(cfs[,2])) +
theme_bw() +
labs(col = "", x = 'standard errors')
```