Assignment 4

Assignment 4 has two purposes:

• To give you more experience with linear regression
• To introduce simulation as a means of studying the behavior of statistical techniques

This assignment is due October 25, 2020 at the end of the day.

Note

In this writeup, I will use capital letters (e.g. $X$) for random variables, and lowercase variables (e.g. $x$) for individual samples (or draws) from those random variables.

You will find the probability notes useful for understanding the derivations in this assignment.

Revision Log

October 19, 2020

Added clarifying notes to Warmup.

Simulation

One common way to understand the behavior of statistical techniques is to use simulation (often called Monte Carlo simulation). In a simulation, we use a psuedorandom number generator to make up data that follows particular patterns (or lack of patterns). We call this data synthetic data.

We then apply a statistical technique, such as correlation coefficient or a linear regression, to the synthetic data, and see how closely its results match the parameters we put in to our simulation. If the analysis reliably estimates the simulation's parameters, we say it recovers the parameters. We can do this many times to estimate that reliability — we can run the simulation 1000 times, for example, and examine the distribution of the error of its parameter estimates to see if it is unbiased, and how broad the errors are.

This technique is commonly used in statistics research (that is, research about statistics itself, rather than research that uses statistics to study other topics) in order to examine the behavior of statistical methods. By simulating samples of different sizes from a population with known parameters, we can compare the results of analyzing those samples with the actual values the statistical method is supposed to estimate. Further, by mapping its behavior over a range of scenarios, we can gain insight into what a statistical technique is likely doing with the particular data we have in front of us.

This is distinct from bootstrapping. In bootstrapping, we are resampling our sample to try to estimate the sampling distribution of a statistic with respect to the population our sample was drawn from; we have actual data, but do not know the actual population parameters. In simulation, we know the population parameters, and do not have any actual data because we make it all up with the random number generator.

Generating Random Numbers

NumPy's Generator class is the starting point for generating random numbers. It has methods for generating numbers from a range of distributions. For more sophisticated distributions, the various distributions in the scipy.stats also support random draws.

Random number generators have a seed that is the starting point for picking numbers. Two identical generators with the same seed will produce the same sequence of values.

We can create a generator with np.random.default_rng:

rng = np.random.default_rng(20201014)

In my class examples, I have been using the current date as my seed. If you do not specify a seed, it will pick a fresh one every time you start the program; for reproducibility, it is advised to pick a seed for any particular analysis. It's also useful to re-run the analysis with a different seed and double-check that none of the conclusions changed.

We can then use the random number generator to generate random numbers from various distributions. It's important to note that random does not mean uniform — then uniform distribution is just one kind of random distribution.

For example, we can draw 100 samples from the standard normal distribution ($\mu =0$, $\sigma =1$):

xs = rng.standard_normal(100)

Warmup: Correlation (10%)

If two variables are independent, their correlation should be zero, right? We can simulate this by drawing two arrays of 100 standard normal variables each, and computing their correlation coefficient:

xs = pd.Series(rng.standard_normal(100))
ys = pd.Series(rng.standard_normal(100))
xs.corr(ys)

Note

This code takes 100 draws (or samples) from the standard normal, twice (once for xs and again for ys).

Mathematically, we write this using $\sim$ as the operator “drawn from”:

$$\begin{array}{rl}x& \sim \mathrm{Normal}(0,1)\\ y& \sim \mathrm{Normal}(0,1)\end{array}$$

Run 1000 iterations of this simulation to compute 1000 correlation coefficients. What is the mean and variance of these simulated coefficients? Plot their distribution. Are the results what you expect for computing correlations of uncorrelated variables?

Tip

What you need to do for this is to run the code example above — that computes the correlation between two 100-item samples — one thousand times. This will draw a total of 200,000 numbers (100 each for x and y, in each simulation iteration).

Repeat the previous simulation, but using 1000 draws per iteration instead of 100. How does this change the mean and variance of the resulting coefficients?

Tip

Now you need to modify the code to draw 1000 normals for x and 1000 normals for y in each iteration. Remember that the example above is drawing 100 normals for each variable.

Remember the covariance of two variables is defined as:

$$\mathrm{Cov}(X,Y)=\mathrm{E}[(X-\mathrm{E}[X])(Y-\mathrm{E}[Y])]=\frac{1}{n}\sum _i(x_i-\overline{x})(y_i-\overline{y})$$

And the correlation is:

$$r=\frac{\mathrm{Cov}(X,Y)}{{\sigma }_X{\sigma }_Y}=\frac{\mathrm{Cov}(X,Y)}{\sqrt{\mathrm{Cov}(X,X)}\sqrt{\mathrm{Cov}(Y,Y)}}$$

If we want to generate correlated variables, we can do so by combining two random variables to form a third:

$$\begin{array}{rl}x& \sim \mathrm{Normal}(0,1)\\ y& \sim \mathrm{Normal}(0,1)\\ z& =x+y\end{array}$$

We can draw them with:

xs = pd.Series(rng.standard_normal(100))
ys = pd.Series(rng.standard_normal(100))
zs = xs + ys

With these variables, we have:

$$\begin{array}{rl}{\sigma }_X& =1\\ \mathrm{E}[X]& =0\\ {\sigma }_Y& =1\\ \mathrm{E}[Y]& =0\\ \mathrm{E}[Z]& =\mathrm{E}[X+Y]=\mathrm{E}[X]+\mathrm{E}[Y]=0\end{array}$$

This last identity is from a property called linearity of expectation. We can now determine the covariance between $X$ and $Z$. As a preliminary, since $X$ and $Y$ are independent, their covariance $\mathrm{Cov}(X,Y)=0$. Further, their independence implies that $\mathrm{E}[XY]=\mathrm{E}[X]\mathrm{E}[Y]$, which from the equations above is $0$.

With that:

$$\begin{array}{rl}\mathrm{Cov}(X,Z)& =\mathrm{E}[(X-\mathrm{E}[X])(Z-\mathrm{E}[Z])]\\ & =\mathrm{E}[X(X+Y)]\\ & =\mathrm{E}[X^2+XY]\\ & =\mathrm{E}[X^2]+\mathrm{E}[XY]\\ & =\mathrm{E}[X^2]\\ & =\mathrm{Var}[X]& \text{since }\mathrm{E}[X]=0\\ & ={\sigma }_X^2=1\end{array}$$

The correlation coefficient depends on $\mathrm{Cov}(X,Z)=1$, $\mathrm{Var}(X)=1$, and $\mathrm{Var}(Z)$. We can derive $\mathrm{Var}(Y)$ as follows:

$$\begin{array}{rl}\mathrm{Var}(Z)& =\mathrm{E}[(Z-\mathrm{E}[Z]{)}^2]\\ & =\mathrm{E}[Z^2]\\ & =\mathrm{E}[(X+Y{)}^2]\\ & =\mathrm{E}[X^2+2XY+Y^2]\\ & =\mathrm{E}[X^2]+2\mathrm{E}[XY]+\mathrm{E}[Y^2]\\ & =\mathrm{Var}(X)+\mathrm{Var}(Y)& \text{since both variables' means are 0}\\ & ={\sigma }_X^2+{\sigma }_Y^2=2\end{array}$$

Therefore we have ${\sigma }_X=1$ (from its distribution), and ${\sigma }_Z=\sqrt{2}$, and so the correlation

$$r_{XZ}=\frac{\mathrm{Cov}(X,Z)}{{\sigma }_X{\sigma }_Z}=\frac{1}{1\cdot \sqrt{2}}=0.707$$

Covariance

You can compute the covariance with the Pandas .cov method. It's instructive to also plot that!

Run 1000 iterations simulating these correlated variables to compute 1000 correlation coefficients (xs.corr(zs)). Compute the mean and variance of these coefficients, and plot their distributions. Does this match what we expect from the analytic results? What happens when we compute correlations of 1000-element arrays in each iteration? What about 10000-element arrays?

Linear Regression (40%)

If we want to simulate a single-variable linear regression:

$$y=\alpha +\beta x+ϵ$$

there are four things we need to control:

• the distribution of $x$
• the intercept $\alpha$
• the slope $\beta$
• the variance of errors ${\sigma }_{ϵ}^2$

Remember that the linear regression model assumes errors are i.i.d. normal, and the OLS model will result in a mean error of 0; thus we have $ϵ\sim \mathrm{Normal}(0,{\sigma }_{ϵ})$. Sampling data for this model involves the following steps:

1. Sample $x$
2. Sample $ϵ$
3. Compute $y=\alpha +\beta x+ϵ$

Let's start with a very simple example: $x$ is drawn from a standard normal, $\alpha =0$, $\beta =1$, and ${\sigma }_{ϵ}^2=1$.

xs = rng.standard_normal(1000)
errs = rng.standard_normal(1000)
ys = 0 + 1 * xs + errs
data = pd.DataFrame({
    'X': xs,
    'Y': ys
})

Fit a linear model to this data, predicting $Y$ with $X$. What are the intercept and slope? What is $R^2$? Are these values what you expect? Plot residuals vs. fitted and a Q-Q plot of residuals to check the model assumptions - do they hold?

Repeat the simulation 1000 times, fitting a linear model each time. Show the mean, variance, and a distribution plot of the intercept, slope, and $R^2$.

Extracting Parameters

The RegressionResults class returned by .fit() contains the model parameters. The .params field has the coefficients (including intercept), and .rsquared has the $R^2$ value:

fit.params['X']

Fit a model to data with $\alpha =1$ and $\beta =4$. Are the resulting model parameters what you expect? How did $R^2$ change, and why? Do the linear model assumptions still hold? What are the distributions of the slope, intercept, and $R^2$ if you do this 1000 times?

Nonlinear Data (15%)

Generate 1000 data points with the following distributions and formula:

$$\begin{array}{rl}x& \sim \mathrm{Normal}(0,1)\\ ϵ& \sim \mathrm{Normal}(0,5)\\ y& =10+5e^x+ϵ\end{array}$$

Fit a linear model predicting $y$ with $x$. How well does the model fit? Do the assumptions hold?

Draw a scatter plot of $x$ and $y$.

Drawing Normals

You can draw from $\mathrm{Normal}(0,5)$ either by using the normal method of Generator, or by drawing an array of standard normals and multiplying it by 5.

Tip

The NumPy function np.exp computes $e^x$.

Repeat with $y=-2+3x^3+ϵ$

Non-Normal Covariates (15%)

Generate 1000 data points with the model:

$$\begin{array}{rl}y& =10+0.1x+ϵ\\ ϵ& \sim \mathrm{Normal}(0,1)\\ x& \sim \mathrm{Gamma}(2,1)\end{array}$$

• Plot the distributions of $X$ and $Y$
• Fit a linear model predicting $y$ with $x$
• How well does this model fit? Do the assumptions hold?

Gamma Distributions

You can draw 1000 from the $Gamma(2,1)$ by:

rng.gamma(2, 1, 1000)

Multiple Regression (10%)

Now we're going to look at regression with two or more independent variables.

We will use the following data generating process:

$$\begin{array}{rl}{x}_1& \sim \mathrm{Normal}(10,2)\\ x_2& \sim \mathrm{Normal}(-2,5)\\ ϵ& \sim \mathrm{Normal}(0,1)\\ y& =1+0.5x_1+3x_2+ϵ\end{array}$$

Tip

To draw from $\mathrm{Normal}(\mu ,\sigma )$, you can draw xs from a standard normal and compute xs * σ + μ.

Fit a linear model y ~ x1 + x2 on 1000 data points drawn from this model. What are the intercept and coefficients from the model? Are they what you expect? Check the model assumptions — do they hold?

Note

You can draw both $x_1$ and $x_2$ simultaneously with:

xs = rng.multivariate_normal([10, -2], [[2, 0], [0, 5]], 100)
# turn into a data frame
xdf = pd.DataFrame(xs, columns=['X1', 'X2'])

The multivariate normal distribution is parameterized by a list (or array) of means, and a positive symmetric covariance matrix defined as follows:

$$\left[\begin{array}{cc}\mathrm{Var}(X_1)& \mathrm{Cov}(X_1,X_2)\\ \mathrm{Cov}(X_2,X_1)& \mathrm{Var}(X_2)\end{array}\right]$$

That is, the diagonals of them matrix are the variances of the individual variables, and the other cells are the covariances between pairs of variables. The example code sets up the following matrix:

$$\left[\begin{array}{cc}2& 0\\ 0& 5\end{array}\right]$$

Correlated Predictors (10%)

Now we're going to see what happens when we have correlated predictor variables. Remember I said those were a problem?

We're going to use the multivariate normal from the hint in the previous part to draw correlated variables $X_1$ and $X_2$ to use as predictors. We will use the following procedure:

1. Draw 1000 samples of variables $X_1$ and $X_2$ from a normal with means $⟨1,3⟩$, variances of 1, and a covariance $\mathrm{Cov}(X_1,X_2)=0.85$:

   xs = rng.multivariate_normal([1, 3], [[1, 0.85], [0.85, 1]], 1000)
   

2. Draw $ϵ\sim \mathrm{Normal}(0,2)$

3. Compute $y=3+2x_1+3x_2+ϵ$

Show a pairplot of our variables $X_1$, $X_2$, and $Y$. What do we see about their distributions and relationships?

Fit a linear regression for y ~ x1 + x2. How well does it fit? Do its assumptions hold?

Run this simulation (drawing 1000 variables and fitting a linear model) 100 times. Show the mean, variance, and distribution plots of the estimated intercepts and coefficients (for x1 and x2).

Repeat the repeated simulation for a variety of different covariances from 0 to 1 (including at least 0, 1, 0.9, and 0.99). Create line plots (or a single plot with multiple colors) that show how the variance of the estimated regression parameters (intercept and $x_1$ and $x_2$ coefficients) change as you increase the correlation (covariance) between $X_1$ and $X_2$.

Reflection (10%)

Write a couple of paragraphs about what you learned from this assignment.

Expected Time

As in A3, I'm providing here some estimates of how long I expect each part might take you.

1. Warmup: 1 hour
2. Linear Regression: 2 hours
3. Nonlinear Data: 30 minutes
4. Non-normal Covariates: 30 minutes
5. Multiple Regression: 1 hour
6. Correlated Predictors: 2 hours
7. Reflection: 1 hour