Simulation of covariates in regression

Question

I am doing a simulation study on a simple linear regression model to see how well OLS estimator performs. I am planning to use coverage probability to assess the estimation of slope parameter using 1000 replications with sample size of 100

Suppose the model is simple as $y_i=0.5+2*x_i+\epsilon_i$

I am not sure which of the following two options to choose in simulation.

1. simulate all 100 $x_i$ values, say from a normal distribution. Then in each of the 1000 replications, use the same 100 $x_i$ values for generation of $y_i$.
2. In each of 1000 replications, first simulate 100 $x_i$ values, and based on these covariate values, generate $y_i$.

Any help is appreciated.

Answer 1

To answer your question in comments. If you were doing a bootstrap instead of a simulation the bootstrap would suggest something like option 1.

If you want to accept the assumption that the xs are fixed and known which is an important assumption for applying OLS then option 1 demonstrates the performance when the assumption of known x is correct. But suppose you are interested in seeing the sensitivity to the assumption of fixed x. Then you want to simulate with random errors generated for both x and y. That would not be option 1 but it would not be option 2 either.

Answer 2

for $i \in \{1, 2, \ldots, 100\}$ let

$$\eta_i = .5 + 2x_i$$

which is the linear predictor and then you can simulate

$$y_i \sim \mathcal{N}(\eta_i, \sigma)$$

If you do $j$ replicates, your question boils down to:

$$y_{ij} \sim \mathcal{N}(\eta_i, \sigma)$$

or

$$y_{ij} \sim \mathcal{N}(\eta_{ij}, \sigma)$$

If you are interested in drawing inference about parameters estimates for a given $\boldsymbol{\eta}$, then I think you want $y_{ij} \sim \mathcal{N}(\eta_i, \sigma)$ where you use the same values for your linear predictor every time (i.e., option 1 in your question).

A simple example

set.seed(10)
x <- rnorm(100, mean = 0, sd = 1)
res1 <- sapply(1:10000, function(i) {
  coef(lm(rnorm(100, mean = .5 + 2 * x, sd = 1) ~ x))
})

set.seed(10)
res2 <- sapply(1:10000, function(i) {
  x <- rnorm(100, mean = 0, sd = 1)
  coef(lm(rnorm(100, mean = .5 + 2 * x, sd = 1) ~ x))
})

rowMeans(res1)
# (Intercept)           x 
#  0.5002331   2.0010302 
apply(res1, 1, sd)
# (Intercept)           x 
#  0.1005374   0.1070732 
rowMeans(res2)
# (Intercept)           x 
#  0.4995853   1.9987682 
apply(res2, 1, sd)
# (Intercept)           x 
#  0.09977722  0.10169706