# Simulation from linear model with additional variables

I want to test the performance of a variable selection method in linear regression with normal errors using simulated data:

$${\bf y}= {\bf X}{\bf \beta} + \epsilon,$$

where, as usual, ${\bf y}$ is $n\times 1$, ${\bf\beta}$ is $p\times 1$, $n>p$, $\epsilon_j\stackrel{ind.}{\sim}N(0,\sigma^2)$, $j=1,\dots,n$.

How can I simulate additional superfluous variables? Is there a benchmark method for adding variables or is it an irrelevant part of the simulation? I was thinking of adding simulated columns, from some arbitrarily chosen distribution, in the design matrix and check whether or not the variable selection method detects the artificial additional variables.

• I have accepted an answer based on the merits of the effort shown. However, the general question about whether or not there is a benchmark method remains (in a formal sense) open. Commented Sep 17, 2014 at 22:04

## 3 Answers

One way would be to simulate all $x_1, x_2, ..., x_p$ together, assign each explanatory variable a coefficient, then simulate the error term $\epsilon$, and finally your dependent variable would just be the sum of the $X'\beta$ and $\epsilon$. Many statistical packages have functions where you can specify the correlation between the $x$ variables, too. In Stata, for instance, that could be achieved with the corr2data command.

Perhaps you are not using Stata but as long as you know the simulation commands in other languages the steps should be the same.

// set a certain number of observations
set obs 1000

// generate the explanatory variables (here we simulate 2 variables from a normal distribution)
gen x1 = rnormal(5,3)
gen x2 = rnormal(9,1)

// generate the error term (here is the simple case where the error is distributed as N(0,1) - for other distributions use the according sampling technique)
gen e = rnormal(0,1)

// generate the dependent variable and assign coefficients to the explanatory variables (0.5 for x1 and 0.9 for x2, for instance)
gen y = 0.5*x1 + 0.9*x2 + e

// run the linear regression of y on x1 and x2
reg y x1 x2


The corr2data command gives you much more options to specify correlations between the variables. So you can see what happens to your model if you have high collinearity between x1 and x2, you can simulate measurment error, correlations with the error, etc. It can also be used to generate a heteroscedastic relationship between one or more of the explanatory variables with the error.

Given the edit of the original question, here is also how to add superfluous variables. For this you would need to specify a correlation matrix before generating the data, for example:

    |       x1       x2       x3        e
----+------------------------------------
x1 |   1.0000
x2 |   0.3000   1.0000
x3 |   0.0100  -0.0000   1.0000
e |   0.0000  -0.0000  -0.0000   1.0000


Which can be achieved via

mat C = (1, 0.3, 0.01, 0 \ 0.3, 1, 0, 0 \ 0.01, 0, 1 , 0 \ 0, 0, 0, 1)
corr2data x1 x2 x3 e, n(1000) means(5 7 13 0) sds(3, 1, 2, 1) corr(C)
corr


where you then make x3 "superfluous" by simply not including it in the construction of y. It's not completely useless because it is correlated with x1, so through the correlation matrix C you can decide how superfluous x3 actually is. Then generating y in the same way as before

// generate the dependent variable
gen y = 0.5*x1 + 0.9*x2 + e

// run the regression with the useless variable
reg y x1 x2 x3


gives the result you wanted.

But this would be the generic set-up to generate such data which should work in any other statistical package in the same way. Presumably other packages have additional/different functions that you can use but the steps done here are a basic way to achieve this.

• Thank you for your answer. If I properly understand, you are explaining how to simulate from the linear model. Now, if I want to add "superfluous" variables (in addition to $x_1,\dots,x_p$, say, $x_{p+1},\dots$), what would you recommend? Commented Sep 17, 2014 at 16:22
• You can simulate other variables and simply add them to the regression. By construction they are then uncorrelated with $y$. The corr2data command, for instance, would also allow generation of "superfluous" variables that are not correlated with $y$ but with the other explanatory variables.
– Andy
Commented Sep 17, 2014 at 16:50

You could also approach to simulating data ($\mathbf{X}$) containing $p$ variables and $n$ observations by using a Cholesky decomposition of an a priori specified $p$-by-$p$ correlation matrix ($\mathbf{R}$) and an $p$-by-$n$ 'noise matrix' ($\mathbf{V}$):

$$\mathbf{X} = \left(\mathbf{R^{*}V}\right)^{\text{T}}$$

where $\mathbf{R^{*}}$ is the Cholesky decomposition of $\mathbf{R}$, and $\mathbf{V} \sim \mathcal{N}\left(0,\sigma\right)$. Your linear relationships would thus be articulated through the correlations defined as elements of $\mathbf{R}$.

References:
Dray, S. (2008). On the number of principal components: A test of dimensionality based on measurements of similarity between matrices. Computational Statistics & Data Analysis, 52(4), 2228–2237.

Linn, R. L. (1968). A Monte Carlo approach to the number of factors problem. Psychometrika, 33(1), 37–71.

1. Generate shocks ($\epsilon_j \sim N(0, \sigma^2))$
2. Let $Y_{truth} = X\beta + \epsilon$
3. Fit $Y_{polluted} = X\beta + P\gamma$
4. Run your selection procedure and hope that $Y_{polluted}$ will reduce to $Y_{truth}$
• Thanks. Is there a standard procedure for choosing $P$? How relevant is this choice? Commented Sep 17, 2014 at 16:34
• In theory the choice of P does not really matter, because you have already generated the ground truth! In this case you know the data generating process so your goal is just to recover the parameters you used for simulating the ground truth. In practice, the choice of P might be important due to numerical stability. For starter, you could choose a P that is as orthogonal as possible to X. Commented Sep 17, 2014 at 16:37
• Thank you. I like your answer, I just wonder if there is a formal argument to justify the irrelevance of the choice of $P$. I think the answer is not straightforward. Commented Sep 17, 2014 at 22:10