R

This simulation can easily be done in R without the use of loops. First, we set the sample size and simulation size of the simulation. Next, we just randomly sample from each of the above distributions using the rnorm, runif, and rexp functions. In each of these functions we sample as many observations we would need across all simulations which is the product of the sample size with the simulation size. We then use the matrix function to take the draws from the random sampling functions and allocate the observations across rows and columns of a matrix. Specifically, we will use the rows as simulations, and columns as samples. Lastly, we just use the apply function for each matrix and apply the mean function to each row (or simulation). From there we can plot the distributions of the sample means.

sample.size <- 50
simulation.size <- 10000

X1 <- matrix(data=rnorm(n=(sample.size*simulation.size), mean=2, sd=5), nrow=simulation.size, ncol=sample.size, byrow=TRUE)
X2 <- matrix(data=runif(n=(sample.size*simulation.size), min=5, max=105), nrow=simulation.size, ncol=sample.size, byrow=TRUE)
X3 <- matrix(data=(rexp(n=(sample.size*simulation.size)) + 3), nrow=simulation.size, ncol=sample.size, byrow=TRUE)

Mean.X1 <- apply(X1,1,mean)
Mean.X2 <- apply(X2,1,mean)
Mean.X3 <- apply(X3,1,mean)

hist(Mean.X1, breaks=50, main='Sample Distribution of Means for Normal Distribution', xlab='Sample Means')

hist(Mean.X2, breaks=50, main='Sample Distribution of Means for Uniform Distribution', xlab='Sample Means')

hist(Mean.X3, breaks=50, main='Sample Distribution of Means for Exponential Distribution', xlab='Sample Means')

The example above shows sampling distributions for the mean when the sample size is 50. Easily changing that piece of the code will allow the user to see these sampling distributions at different sample sizes.

As you can see above, the sampling distribution for the sample means is approximately normally distributed, regardless of the population distribution, as long as the sample size is big enough.

SAS

This simulation can easily be done in SAS with a couple of loops in a DATA step. First, we set the sample size and simulation size of the simulation with MACRO variables using the % symbol. Next, we open two DO loops to loop over the simulation and sample size. Then we just randomly sample from each of the above distributions using the RAND function with the Normal, Uniform, and Exponential options. We then use the OUTPUT statement to make sure that the results from each iteration in the loop are stored as a separate row in the dataset. Lastly, we just use the PROC MEANS procedure to calculate the means across each variable (using the VAR statement) for each simulation (using the BY statement). We use the OUTPUT statement and OUT option to store the resulting means into separate datasets as different columns - one column for each simulation. From there we can plot the distributions of the sample means using three separate SGPLOT procedure calls.

%let Sample_Size = 50;
%let Simulation_Size = 10000;

data CLT;
    do sim = 1 to &Simulation_Size;
        do obs = 1 to &Sample_Size;
            call streaminit(12345);
            X1 = RAND('Normal', 2, 5);
            X2 = 5 + 100*RAND('Uniform');
            X3 = 3 + RAND('Exponential');
            
            output;
        end;
    end;
run;

proc means data=CLT noprint mean;
    var X1 X2 X3;
    by sim;
    output out=Means mean(X1 X2 X3)=Mean_X1 Mean_X2 Mean_X3;
run;

proc sgplot data=Means;
    title 'Sampling Distribution of Means';
    title2 'Normal Distribution';
    histogram Mean_X1;
    keylegend / location=inside position=topright;
run;

proc sgplot data=Means;
    title 'Sampling Distribution of Means';
    title2 'Uniform Distribution';
    histogram Mean_X2;
    keylegend / location=inside position=topright;
run;

proc sgplot data=Means;
    title 'Sampling Distribution of Means';
    title2 'Exponential Distribution';
    histogram Mean_X3;
    keylegend / location=inside position=topright;
run;

The example above shows sampling distributions for the mean when the sample size is 50. Easily changing that piece of the code will allow the user to see these sampling distributions at different sample sizes.

As you can see above, the sampling distribution for the sample means is approximately normally distributed, regardless of the population distribution, as long as the sample size is big enough.

R

This simulation is a little more drawn out. First, we set the parameters of interest - sample size, simulation size, coefficients for the regression model, and the correlation between the X’s. Next, we calculate the Cholesky decomposition of our correlation matrix to help with the correlating of data. From there we simulate three variables, \(X_1\), \(X_2\), and \(X_{2u}\). The first \(X_2\) we will correlate with \(X_1\), but we will leave \(X_{2u}\) uncorrelated with \(X_1\). Lastly, we correlate \(X_1\) and \(X_2\) as described in the above sections on correlated variables.

sample.size <- 50
simulation.size <- 10000
beta0 <- -13
beta1 <- 1.21
beta2 <- 3.45
correlation_coef <- 0.6

R <- matrix(data=cbind(1,correlation_coef, correlation_coef, 1), nrow=2)
U <- t(chol(R))

X1 <- rnorm(sample.size*simulation.size, mean = 0, sd = 1)
X2 <- rnorm(sample.size*simulation.size, mean = 0, sd = 1)
X2u <- rnorm(sample.size*simulation.size, mean = 0, sd = 1)

Both <- cbind(standardize(X1), standardize(X2))
X1X2 <- U %*% t(Both)
X1X2 <- t(X1X2)

Now we need to create our models. In other words, we need to simulate our target variables under the conditions described above. We have two models to simulate. The first is Y.c where we have \(Y\) as a function of \(X_1\), \(X_2\), and our error, \(\varepsilon\), simulated from a normal distribution using rnorm. The second is Y.u where we have \(Y\) as a function of \(X_1\), \(X_2u\), and our error, \(\varepsilon\), simulated from a normal distribution using rnorm.

Next, we just create columns for each simulation of our variables in datasets to loop across.

# Time to create the two target variables - one with correlated X's, one without. #
all_data <- data.frame(X1X2)

all_data$Y.c <- beta0 + beta1*all_data$X1 + beta2*all_data$X2 + rnorm(sample.size*simulation.size, mean = 0, sd = 1.5)
all_data$Y.u <- beta0 + beta1*all_data$X1 + beta2*X2u +rnorm(sample.size*simulation.size, mean = 0, sd = 1.5)

# Getting all the data in separate columns to loop across. #
X1.c <- data.frame(matrix(all_data$X1, ncol = simulation.size))
X2.c <- data.frame(matrix(all_data$X2, ncol = simulation.size))
X2.uc <- data.frame(matrix(X2u, ncol = simulation.size))
Y.c <- data.frame(matrix(all_data$Y.c, ncol = simulation.size))
Y.uc <- data.frame(matrix(all_data$Y.u, ncol = simulation.size))

Now we get to run the regressions. We use lapply to apply the lm function (for building linear regressions) across all three types of models. The all_results_c object has the linear regressions leaving out a correlated \(X_2\). The all_results_uc object has the linear regressions leaving out a uncorrelated \(X_2u\). The all_results_ols object has the correct linear regressions leaving out no variables. From there, we use the sapply function to apply the coef function across each of our linear regressions to get the coefficients from the models. Lastly, we just plot the densities of the coefficients from each of the three scenarios using the density function.

# Running regressions across every column. #
all_results_c <- lapply(1:simulation.size, function(x) lm(Y.c[,x] ~ X1.c[,x]))
all_results_uc <- lapply(1:simulation.size, function(x) lm(Y.uc[,x] ~ X1.c[,x]))
all_results_ols <- lapply(1:simulation.size, function(x) lm(Y.c[,x] ~ X1.c[,x] + X2.c[,x]))

# Gathering and plotting the distribution of the coefficients. #
betas_c <- sapply(all_results_c, coef)
betas_uc <- sapply(all_results_uc, coef)
betas_ols <- sapply(all_results_ols, coef)

par(mfrow=c(3,1))
plot(density(betas_ols[2,]), xlim = c(-1,5), main = "Slope Histogram - OLS with Correct Model")
abline(v = 1.21, col="red", lwd=2)

plot(density(betas_c[2,]), xlim = c(-1,5), main = "Slope Histogram - Correlated X2 Not in Model")
abline(v = 1.21, col="red", lwd=2)

plot(density(betas_uc[2,]), xlim = c(-1,5), main = "Slope Histogram - Uncorrelated X2 Not in Model")
abline(v = 1.21, col="red", lwd=2)

par(mfrow=c(1,1))

From the plots above we can see the impact of leaving out important variables that should be in our model. If our model leaves out a predictor variable that is correlated with another predictor variable, the coefficient of the predictor variable in the model is biased. Here we can see, the vertical lines showing the true value of the coefficient. When the omitted variable is correlated, our distribution of the coefficient of \(X_1\) is not centered around its true value. Even if the omitted variable isn’t correlated, we can see that the standard error (the spread) of the distribution of the coefficient for the variable \(X_1\) is still wider than it would be under the circumstance where we specified the model correctly.

Another thing we can check is how many times we incorrectly, fail to reject our null hypothesis that \(X_1\) isn’t a significant variable in our model. We know it is an important variable since we created the model. However, let’s look at the p-values for the coefficient on \(X_1\) for each of the simulations. We just use the sapply function to apply the summary function and extract the coefficients element from each linear regression. Then we just sum up how many times our p-value was above 0.05. This should NOT be the case too often since we know \(X_1\) should be in our model.

pval_c <- sapply(all_results_c, function(x) summary(x)$coefficients[2,4])
pval_uc <- sapply(all_results_uc, function(x) summary(x)$coefficients[2,4])
pval_ols <- sapply(all_results_ols, function(x) summary(x)$coefficients[2,4])

(sum(pval_ols > 0.05)/10000)*100

[1] 1.36

(sum(pval_c > 0.05)/10000)*100

[1] 0

(sum(pval_uc > 0.05)/10000)*100

[1] 40.9

As you can see, in the correctly specified model, we incorrectly say \(X_1\) is not statistically significant only 1.39% of the time. However, when we left out our \(X_2\) variable and it wasn’t correlated, we make this mistake 40.84% of the time.

Omitting important variables is a problem. If those variables are correlated with other inputs that are in the model, then our inputs in the model have biased estimations. However, even if they are not correlated and therefore unbiased, we make many more mistakes about our model input’s significance in the model.

SAS

This simulation is a little more drawn out. First, we set the parameters of interest - sample size, simulation size, coefficients for the regression model, and the correlation between the X’s using MACRO variables. Next, we build in the correlation structure between \(X_2\) and \(X_2u\) using the same approach as descibed in the correlation section above using DATA steps and PROC MODEL.

%let Sample_Size = 50;
%let Simulation_Size = 10000;
%let Beta0 = -13;
%let Beta1 = 1.21;
%let Beta2 = 3.45;
%let Correlation_Coef = 0.6;

   /* Creating correlation matrices */
data Corr_Matrix;
    do i = 1 to &Simulation_Size;
        _type_ = "corr";
        _name_ = "X_1";

        X_1 = 1.0; 
        X_2 = &Correlation_Coef;

        output;
        _name_ = "X_2";

        X_1 = &Correlation_Coef; 
        X_2 = 1.0;

        output;
    end;
run;

   /* Creating initial values as place holders */
data Naive;
    do i = 1 to &Simulation_Size;
        X_1 = 0;
        X_2 = 0;
        output;
    end;
run;

   /* Generate the correlated X_1 and X_2 */
proc model noprint;
  /* Distributiuon of X_1 */
  X_1 = 0;
  errormodel X_1 ~ Normal(1); 

  /*Distribution of X_2*/
  X_2 = 0;
  errormodel X_2 ~ Normal(1); 

  /* Generate the data and store them in the dataset work.Correlated_X */
  solve X_1 X_2 / random = &Sample_Size sdata = Corr_Matrix
  data = Naive out = Correlated_X(keep = X_1 X_2 i _rep_);
    by i;

  run;
quit;

data Correlated_X;
    set Correlated_X;
    by i;
    if first.i then delete;
    rename i = sim;
    rename _rep_ = obs;
run;

Now we need to create our models. First, we will simulate our errors from the model, \(\varepsilon\) and our uncorrelated \(X_2u\) variable. Then we can merge all of our simulated data together in a DATA step with the MERGE statement. With this new combined dataset we have columns of all of our simulated variables across many simulated rows. All we then need to do is combine the columns into our target variables \(Y\) and \(Yu\) for the model with correlated and uncorrelated inputs respectively.

Next, we can easily build .

/* Generate Linear Regression Model Data */
data simulation;
do sim = 1 to &Simulation_Size;
    do obs = 1 to &Sample_Size;
        call streaminit(112358);
        err = sqrt(2.25)*rand('normal');
        X_2u = sqrt(1)*rand('normal');

        output;
    end;
end;
run;

data all_data;
    merge simulation Correlated_X;
    by sim obs;

    Y = &Beta0 + &Beta1*X_1 + &Beta2*X_2 + err;
    Yu = &Beta0 + &Beta1*X_1 + &Beta2*X_2u + err;
run;

Now we get to run the regressions. We use PROC REG to build all of our models by the variable sim using the BY statement so we get results for each simulation’s rows. We store the output from all these regressions into a dataset called all_results using the OUTEST = option. We are complete with the simulation. Next we just clean up the results with some DATA step manipulation for labeling as well as a PROC SORT to make the visuals displayed in PROC SGPANEL nicer to see. The options in the plotting procedure aren’t covered in detail here as that is not the focus.

/* Run all regressions and store the output */
proc reg data = all_data outest = all_results tableout noprint;
    OLSC: model Y = X_1 X_2;
    Corr: model Y = X_1;
    UnCorr: model Yu = X_1;
    by sim;
run;

/* Give meaningful values to the _model_ variable */
data all_results;
    set all_results;
    label _model_="Model";
    if upcase(_model_) eq "Corr" then _model_="Omit Correl. X";
    else if upcase(_model_) eq "UnCorr" then _model_="Omit Uncorrel. X"; 
    else if upcase(_model_) eq "OLSC" then _model_="Correct Model";
run;

/* Normal plots and univariate statistics for each variable/model */
proc sort data=all_results;
    by _model_;
run;

/* Create a panel of plots, by model, to contrast normal vs. empirical distr. AKA Make things look pretty... */
/* Also add a reference line that highlights the true parameter value */
proc sgpanel data = all_results(where=(_type_ eq 'PARMS'));
    panelby _model_ / layout=ROWLATTICE uniscale=column columns=1 onepanel;
    density X_1 / type=kernel;
    refline &Beta1 / axis=x label="&Beta1" labelattrs=(weight=bold);
run;

From the plots above we can see the impact of leaving out important variables that should be in our model. If our model leaves out a predictor variable that is correlated with another predictor variable, the coefficient of the predictor variable in the model is biased. Here we can see, the vertical lines showing the true value of the coefficient. When the omitted variable is correlated, our distribution of the coefficient of \(X_1\) is not centered around its true value. Even if the omitted variable isn’t correlated, we can see that the standard error (the spread) of the distribution of the coefficient for the variable \(X_1\) is still wider than it would be under the circumstance where we specified the model correctly.

Omitting important variables is a problem. If those variables are correlated with other inputs that are in the model, then our inputs in the model have biased estimations. However, even if they are not correlated and therefore unbiased, we still have wider standard errors (spread) in out estimates which would lead to more incorrect conclusions in hypothesis tests for the included input variable.

R

Similar to the omitted variable bias example above, we will need to generate predictor variables first and then calculate our target variable from the two that should be in the model. Our first step is to generate a matrix of 100 observations (rows) of 8 variables (columns) all with normal distributions with mean 0 and standard deviation of 1. We will save this as a data frame Fake. Next, we have the Err object that is our model error term, \(\varepsilon\), that has a normal distribution with mean of 0 and standard deviation of 8. Now we create our target variable \(Y\) as defined above and combine all these columns into one data frame Fake.

set.seed(1234)
Fake <- data.frame(matrix(rnorm(n=(100*8)), nrow=100, ncol=8))
Err <- rnorm(n=100, mean=0, sd=8)
Y <- 5 + 2*Fake$X2 - 3*Fake$X8 + Err
Fake <- cbind(Fake, Err, Y)

The next step is to shuffle the target variable. We will create another matrix with 100 observations (rows) and 1000 columns - one column for each shuffle. With a for loop we can randomly sample from a uniform distribution using the runif function. We will then order our Y column by these uniform random numbers. This will essentially randomly order the target variable. We do this repeatedly and store these in each column of our 1000 shuffle columns. We rename these columns (Y.1, Y.2, …) and combine them all into one new data frame.

sim <- 1000

Y.Shuffle <- matrix(0, nrow=100, ncol=sim)
for(j in 1:sim){
  Uniform <- runif(100)
  Y.Shuffle[,j] <- Y[order(Uniform)]
}

Y.Shuffle <- data.frame(Y.Shuffle)
colnames(Y.Shuffle) <- paste('Y.',seq(1:sim),sep="")

Fake <- data.frame(Fake, Y.Shuffle)

Now we just run our 1001 and linear regressions. The for loop will run a linear regression using the lm function on all the shuffled columns. The summary function’s adj.r.squared element is then stored for each target shuffled model. Lastly, we also run the true regression model and get its adjusted \(R^2\) using the same procedure. This makes our 1001 adjusted \(R^2\) metrics. Let’s visualize these adjusted \(R^2\) values on a histogram.

R.sq.A <- rep(0,sim)
for(i in 1:sim){
  R.sq.A[i] <- summary(lm(Fake[,10+i] ~ Fake$X1 + Fake$X2 + Fake$X3 + Fake$X4
                          + Fake$X5 + Fake$X6 + Fake$X7 + Fake$X8))$adj.r.squared
}
True.Rsq.A <- summary(lm(Fake$Y ~ Fake$X1 + Fake$X2 + Fake$X3 + Fake$X4
                         + Fake$X5 + Fake$X6 + Fake$X7 + Fake$X8))$adj.r.squared

hist(c(R.sq.A,True.Rsq.A), breaks=50, col = "blue", main='Distribution of Adjusted R-Squared Values', xlab='Adjusted R-Squared')
abline(v = True.Rsq.A, col="red", lwd=2)
mtext("True Model", at=True.Rsq.A, col="red")

Our model out performs all of the other models. However, despite not expecting any good results from the shuffled data, we see some adjusted \(R^2\) values that almost reach 0.2. Again, this is due to random chance. This allows us to test our model to see if we just got a data set that is lucky, or if we really have some interesting findings.

Let’s see how many variables were significant in each of our shuffled models. Remember, the data has been completely shuffled so there shouldn’t be any significant variables. TO do this, use a similar for loop as above, but instead of looking at the adj.r.squared element of the summary function, we look at the \(4^{th}\) column of our coefficient element - the p-values for our variables. We then count the number of times the p-value was below a significance level of 0.05.

P.Values <- NULL
for(i in 1:sim){
  P.V <- summary(lm(Fake[,10+i] ~ Fake$X1 + Fake$X2 + Fake$X3 + Fake$X4
                    +                           + Fake$X5 + Fake$X6 + Fake$X7 + Fake$X8))$coefficients[,4]
  P.Values <- rbind(P.Values, P.V)
}

Sig <- P.Values < 0.05

colSums(Sig)

(Intercept)     Fake$X1     Fake$X2     Fake$X3     Fake$X4     Fake$X5 
       1000          51          46          53          53          43 
    Fake$X6     Fake$X7     Fake$X8 
         53          46          50 

table(rowSums(Sig)-1)

  0   1   2   3   4 
675 264  53   7   1 

hist(rowSums(Sig)-1, breaks=25, col = "blue", main='Count of Significant Variables Per Model', xlab='Number of Sig. Variables')

Most of the time we see that no variables were significant, but out of the 1000 shuffles, each variable was significant about 50 times. This isn’t surprising as the goal of a significance level of 0.05 in a hypothesis test is to make an error 5% of the time. However, sometimes these errors happened together. Seven out of the 1000 shuffled models had three significant variables in the model. One even had four significant variables. This is all due to random chance.

SAS

Similar to the omitted variable bias example above, we will need to generate predictor variables first and then calculate our target variable from the two that should be in the model. Our first step is to generate a dataset of 100 observations (rows) of 8 variables (columns) all with normal distributions with mean 0 and standard deviation of 1. Next, we have the Err object that is our model error term, \(\varepsilon\), that has a normal distribution with mean of 0 and standard deviation of 8. Now we create our target variable \(Y\) as defined above and combine all these columns into one dataset Fake using the DATA step and the RAND function with the Normal option. Don’t forget the OUTPUT statement to store your DO loop results at each iteration.

data Fake;
    do obs = 1 to 100;
        call streaminit(12345);
        X1 = RAND('Normal', 0, 1);
        X2 = RAND('Normal', 0, 1);
        X3 = RAND('Normal', 0, 1);
        X4 = RAND('Normal', 0, 1);
        X5 = RAND('Normal', 0, 1);
        X6 = RAND('Normal', 0, 1);
        X7 = RAND('Normal', 0, 1);
        X8 = RAND('Normal', 0, 1);
        Err = RAND('Normal', 0, 8);

        Y = 5 + 2*X2 - 3*X8 + Err;
            
        output;
    end;
run;

The next step is to shuffle the target variable. We will do this using a MACRO we create called shuffle. This MACRO will use a DATA step to create a random sample from a uniform distribution using the RAND function with the Uniform option. We will then order our Y column by these uniform random numbers using PROC SORT. This will essentially randomly order the target variable. We use another DATA step to merge this new randomly shuffled dataset with our original one using the MERGE statement. We then use PROC DATASETS to delete this random shuffle dataset to not clog up our memory. Next we create another MACRO called sim_shuffle that runs the shuffle MACRO a specified number of times. We will run ours 500 times.

%macro shuffle(sim);
    data shuffle_∼
        set Fake;
        Uni = RAND('Uniform');
        Y_&sim = Y;
        keep Y_&sim Uni;
    run;

    proc sort data=shuffle_∼
        by Uni;
    run;

    data Fake;
        merge Fake shuffle_∼
        drop Uni;
    run;

    proc datasets noprint;
        delete shuffle_∼
    run;
    quit;
%mend shuffle;

%macro sim_shuffle(sim);
    %do i = 1 %to ∼
        %shuffle(&i);
    %end;
%mend sim_shuffle;

%sim_shuffle(500);

Now we just run our 1001 and linear regressions. The PROC REG procedure will run a linear regression on the original target \(Y\) and on all the shuffled columns. We will use the SELECTION = BACKWARD option to backward select our significant variables with a significance level of 0.05. We store the adjusted \(R^2\) values using the OUTEST option. Let’s visualize these adjusted \(R^2\) values on a histogram with PROC SGPLOT.

proc reg data=Fake noprint outest=Rsq_A;
    model Y Y_1-Y_500 = X1-X8 / selection=backward slstay=0.05 adjrsq;
run;
quit;

proc sgplot data=Rsq_A;
    title 'Distribution of Adjusted R-Squared';
    histogram _ADJRSQ_;
    refline 0.178 / axis=x label='True Model' lineattrs=(color=red thickness=2);
    keylegend / location=inside position=topright;
run;

Our model out performs all of the other models. However, despite not expecting any good results from the shuffled data, we see some adjusted \(R^2\) values that almost reach 0.2. Again, this is due to random chance. This allows us to test our model to see if we just got a data set that is lucky, or if we really have some interesting findings.

Let’s see how many variables were significant in each of our shuffled models. Remember, the data has been completely shuffled so there shouldn’t be any significant variables. TO do this, use a similar for loop as above, but instead of looking at the adj.r.squared element of the summary function, we look at the \(4^{th}\) column of our coefficient element - the p-values for our variables. We then count the number of times the p-value was below a significance level of 0.05.

proc sgplot data=Rsq_A;
    title 'Distribution of Number of Significant Variables';
    histogram _IN_;
run;

Most of the time we see that no variables were significant, but out of the 1000 shuffles, each variable was significant about 50 times. This isn’t surprising as the goal of a significance level of 0.05 in a hypothesis test is to make an error 5% of the time. However, sometimes these errors happened together. Seven out of the 1000 shuffled models had three significant variables in the model. This is all due to random chance.