STATS506-Problem Set 3 Solved

Question 1 [60 points]
First, repeat question 3 parts a-c from problem set 1 using data.table for all computations and data manipulations.

Then, formulate and state a question answerable using the RECS data. Your question should be similar in scope to (one of) parts a-c above and should rely on one or more variables not previously used. Answer your question (using data.table) and provide supporting evidence in the form of nicely formatted graphs and/or tables.

Question 2 [25 points]
In this question you will design a Monte Carlo study in R to compare the performance of different methods that adjust for multiple comparisons (https://xkcd.com/882/). You can read more about each of these methods by referring to help(p.adjust) in R and the references listed there.

Throughout this question, let n= 1000, p = 100 and

                                                                                       βi = { 01        elsei ∈ .{1, … , 10},

Let X ∈ℝn×p with X ∼ N(0p, Σ) and Y ∼ N(Xβ, σ2Ip) where I is a p × p diagonal matrix.

a.   Write a function that accepts matrices X and beta and returns a p by mc_rep matrix of p-values corresponding to p-values for the hypothesis tests:

H0 : βi = 0, H1 : βi ≠ 0.

In addition to X and beta your function should have arguments sigma ( σ) and mc_rep controlling the error variance of

Y and number of Monte Carlo replicates, respectively. Your function should solve the least-squares problems using the QR decomposition of

X′X. This decomposition should only be computed once each time your function is called. i. Refer to the course notes to find β̂.

ii.    Use Y and Ŷ = Xβ̂ to estimate the error variance for each Monte Carlo trial m:

n

                                                                                                     σ2m̂ =    n−1 p ∑i= 1 (Yim − Ŷim)2

iii.   Use the result from ii and the QR decomposition to find the variance of β̂i, vi = σ2̂ (X′X)− 1.

[Note: you will need to do some algebra to determine how to compute (X′X)− 1 using Q and R.

Or you can use the function chol2inv() .] iv. Form Zi = β̂i/vi and find p = 2(1 − Φ− 1(|Zi|)).

Test your function with a specific X and Y by comparing to the output from appropriate methods applied to the object returned by lm(Y ~ 0 + X) . It’s okay if there is some finite precision error less than ~1e-3 in magnitude. Hint: use set.seed() to generate the same Y inside and outside the scope of the function for the purpose of testing.

b.   Choose Σ and σ as you like. Use the Cholesky factorization of Σ to generate a single X. Pass X, β, and σ to your function from the previous part.

c.    Write a function evaluate that takes a set of indices where β≠ 0 and returns Monte Carlo estimates for the following quantities:

 The family wise error rate

The false discovery rate

The sensitivity

The specificity

See this page (https://en.wikipedia.org/wiki/Sensitivity_and_specificity#Sensitivity_index) for additional details.

d.   Apply your function from the previous part to the matrix of uncorrected P-values generated in part B.

Use the function p.adjust() to correct these p-values for multiple comparisons using ‘Bonferroni’, ‘Holm’, ‘BH’ (Benjamini-Hochberg), and ‘BY’. Use your evaluate() function for each set of adjusted p-values.

e.   Produce one or more nicely formatted graphs or tables reporting your results. Briefly discuss what you found.

Question 3 (Optional) [25 points, 5 points each]
This is a bonus question related to problem 6 from the midterm. First, review the script written in Stata available here (https://github.com/jbhender/Stats506_F18/tree/master/solutions/PS3). In this question, you will work through various options for translating this analysis into R. You may submit all or some of these, but each part must be entirely correct to earn the points listed.

a.   Write a translation using data.table for the computations. [5 pts]

b.   Write a function to compute the univariate regression coefficient by group for an arbitrary dependent, independent, and grouping variables. Use data.table for computations within your function. Test your function by showing it produces the same results as in part a. [10 pts]

c.    Compute the regression coefficients using the dplyr verb summarize_at() . [5 pts]

d.   Write a function similar to the one in part b to compute arbitrary univariate regression coefficients by group. Use dplyr for computations within your function. You should read the “Programming with dplyr” vignette at vignette('programming', 'dplyr') before attempting this. Warning: this may be difficult to debug! [10 points]