2
$\begingroup$

I am doing a "Statistical learning in R" online course which requires me to do exercises in R and right now, i'm stuck at a part that requires me to use the bootstrap.

I'm not sure whether it is allowed to ask about homework problems on here, but i have been looking around different forums and tutorials for a while and i can't really put my head around this.

So the original question is

Now, use the (standard) bootstrap to estimate s.e.($ \hat{\beta_1} $) . To within 10%, what do you get?

A dataset is given and in previous questions i was required to estimate the standard error of $\beta_1$ from a linear model on the data. The set Xy consists of three variables with 1000 entries each (X1, X2 and Y).

In my attempt of a solution, i have tried to stay close to the original code used in the Video lesson.

library(ISLR) 
library(boot)
alpha.fn=function(data,index){
X=data$X[index]
Y=data$Y[index]
return((var(Y)-cov(X,Y))/(var(X)+var(Y)-2*cov(X,Y)))}
#"renaming my variables to fit the alpha function"
Xy$X=Xy$X1;
Xy$Y=Xy$y;
set.seed(1)
alpha.fn(Xy,sample(100,100,replace=T))
boot(Xy,alpha.fn,R=1000)

Now i don't know whether this is the "standard" bootstrap i intend to use or how to pass the function that i want to estimate the std. error to within 10%. does it mean i just want it to 0.1 accuracy?

$\endgroup$
  • $\begingroup$ Is this using the ISLR textbook? If so, and you know the exercise number, there are online solutions here $\endgroup$ – MichaelChirico Sep 20 '16 at 15:23
  • 1
    $\begingroup$ If it is about learning to bootstrap, I'd compute the statistic by lm and then focus on bootstrap without using boot. It is not difficult, just a for loop. No idea what is meant by "within 10%", it is strange. $\endgroup$ – Michael M Sep 20 '16 at 17:46
  • 1
    $\begingroup$ I don't think you are using the right formula. In the video, the alpha formula is only used for the investment allocation example. But in this data $X_y$, the formula is linear regression ($y=\beta_0 +\beta_1X_1 +\beta_2X_2$). $\endgroup$ – Tung Nguyen Jun 3 '17 at 11:02
  • 1
    $\begingroup$ @MichaelM - It's how much you can round off the s.e., for example, if the estimate is 3.2, you could report 3, but if it's 3.4, you'd want to report 3.5 (3 is not within 10% of 3.4). If the OP just reports two significant digits, that would do. $\endgroup$ – jbowman Jun 3 '17 at 14:52
  • $\begingroup$ @jbowman: oh, that makes sense. Thx for clarification $\endgroup$ – Michael M Jun 3 '17 at 22:59
1
$\begingroup$

First of all, Tung Nguyen is correct; you have to estimate the full model and extract the relevant estimate (for $\beta_1$), not just run a subset of the model using only $X_1$ on the r.h.s.

Following up on Michael M's suggestion, which is also how I would approach it, we get the following code, with some of R's printout removed for brevity:

# Run basic regression, calculate coefs and s.e.s
summary(lm(y~x1+x2, data=df))

* blah blah blah *

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.03242    0.03213   32.14   <2e-16 ***
x1           0.97044    0.03214   30.19   <2e-16 ***
x2           1.01010    0.03291   30.69   <2e-16 ***
---

# Bootstrap (1000 replications)
beta_hat <- rep(0,1000)
for (i in 1:length(beta_hat)) {
    local_df <- df[sample(1000, replace=TRUE),]
    beta_hat[i] <- coef(lm(y~x1+x2, data=local_df))["x1"]
}
sd(beta_hat)
[1] 0.03258774

Our bootstrap estimate of the standard error is, as one would expect given the large bootstrap replication count and Normally-distributed errors, very close to the calculation made by lm.

You could shorten the R code a little by rewriting the two statements inside the loop as one"

beta_hat[i] <- coef(lm(y~x1+x2, data=df[sample(1000, replace=TRUE),]))["x1"]

but this is likely less clear for expository purposes!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.