# Why doesn't clustering bias the standard errors in this very simple simulation?

I simulate a very simple model with one regressor and a random intercept: $$y_{ij} = \beta_0 + \beta_1x_{1,ij} + u_j + \epsilon_{ij}$$ I run model <- lm(y ~ x1) a thousand times and store each time the estimated $$\hat{\beta_1}$$ and the estimated standard error of that slope:

myBeta <- model$coefficients[2] mySE <- coef(summary(model))[2,2]  As expected, the higher the variance of the random intercept $$u_j$$, the higher the standard deviation of my vector of 1,000 estimated slopes (i.e. sd(myBetas)). But, comparing sd(myBetas) against mean(mySEs), I would have expected the latter to be smaller and the difference between the two to get worse with higher variances of $$u_j$$. Instead, they are basically the same number (i.e the estimated standard errors are not biased). Isn't clustering supposed to cause underestimation of the standard errors? Here is the main part of the code with $$\sigma_{u_j} = 3$$:  group <- rep(1:10, each=100) # Cluster indicator for 10 clusters of 100 units x1 <- rnorm(1000) beta.0 <- 1 beta.1 <- 1.5 y <- beta.0 + beta.1*x1 + rnorm(10,0,3)[group] + rnorm(1000) model <- lm(y ~ x1) myBeta <- model$coefficients[2]
mySE <- coef(summary(model))[2,2]

• Because it is simulation, so should have no business secret. Could you put the sample size, number of cluster, values of parameters into your question? – user158565 Oct 31 '18 at 14:09
• @a_statistician done! – Alvaro Fuentes Oct 31 '18 at 14:31
• I verified your conclusion in SAS: after 1000 replicates, mean of estimated variance of$\hat \beta_1$ = 0.0091578 with SD = 0.0041452, variance of estimated $\beta_1$ = 0.0092646. – user158565 Oct 31 '18 at 15:09
• When I added random intercept in the model, I got variance of estimate = 0.000999939 and mean of estimated variance = 0.0010122 with variance = 4.3366782E-9 – user158565 Oct 31 '18 at 22:34