Question 3: Multiple Hypothesis Testing
Problems happen when p-values are used to determine statistical significance for several covariates at the same time. This is because the guarantee on Type I error only applies to each covariate individually, but people often forget this when they report statistical results. In statistical jargon, this is called “failing to control for multiple hypothesis testing”.
As one example of this, imagine having a dataset that contains DNA sequences and a continuous health outcome for many people. DNA sequences are being collected by companies like 23andme that sell ancestry services and, while many health outcomes data are protected by privacy laws in many countries, it isn’t hard to imagine DNA sequences and health outcomes data being brought together at scale. Some health outcomes, such as life span, are publicly available for some individuals on the web in news reports, social media, and as a part of the public record in some states.
Suppose that a person wishes to use the DNA sequence to predict a continuous health outcome. We will represent the DNA sequence with a large number of covariates, each taking the values 0 and 1. (For example, each binary variable can represent the presence of a particular genetic sequence at a particular place in a particular gene.)  Each individual in the dataset is represented by a number i and each covariate a number j.
In practice, unless the covariates were chosen extremely carefully, most covariates would have no predictive power for the health outcome. In this context, basing a judgment of statistical significance on standard p-values (the kind we have learned about in class and that are reported by statsmodels) can lead a person to incorrectly believe that many covariates do have predictive power.
To generate a dataset, use the following code. Note that none of the predictors in the matrix x have an impact on the outcome y.
n = 1000 # number of datapoints
	p = 100 # number of predictors
	# x is a random nxp matrix of 0s and 1s
	# Each entry is 1 with probability 0.1
	x = 1.0*(np.random.rand(n,p) <= 0.1) 
	y = 15 + np.random.randn(n)
(a) Using the method for generating data above, how many j (from j=0 to j=p) have βj different from 0?  Hint: In this class, we use the notation βj to refer to numbers used to generate data (either simulated data or by nature to generate real data) and we use the notation j to refer to estimates produced by linear regression.  βj never refers to estimates. To answer this question, which is about βj, you should not be looking at the output of linear regression.
(b) Using a for loop, generate 1000 independent random datasets using the code above. For each dataset:
A. fit linear regression, assuming an intercept β0 and the p=100 covariates in x
B. check whether the p-value for β1 is less than 0.05.
C. count the number of p-values that are less than 0.05. This will be a number between 0 and 101, since there are 101 p-values returned when fitting linear regression to this dataset. This is the number of covariates that we would believe to be predictive of the health outcome if we used the p-values in a naive way. 
Note for items B & C: You can get the p-values as an array from a fitted linear model using model.pvalues.
What fraction of your 1000 datasets resulted in the p-value for β1 being less than 0.05? You should see that this is close to 0.05. This is because p-values work to control the Type I error for an individual parameter.
Make a histogram of item C computed above: the number of p-values less than 0.05 in each of the 1000 datasets.
You should see that the number of p-values less than 0.05 is often larger than the number you answered in part (a). This is because we have 100 chances to incorrectly say that a variable is statistically significant, and so the number of times we get at least one of these 100 incorrect is much larger than 0.05