boot() generating an error on replacement - R

Question

I've written a couple of functions for retrieving statistics (coefficients and p-values) from an lm object, to be bootstrapped upon. The coefficient one works; the p-value one is failing with error:

Error in boot(data = data, statistic = bs_p, R = 1000) : 
  number of items to replace is not a multiple of replacement length

I now believe the error is related to the inclusion of a factor variable. Attempting to recreate the problem with easily reproducible data.

    L3 <- LETTERS[1:3]
    data <- data.frame(cbind(x = 20:69, y = 1:50), fac = sample(L3, 50, replace = TRUE))
    bs_p <- function (data, i) {
      d <- data[i,]
      fit <- lm (d$y~d$x*d$fac, data=d)
      return(summary(fit)$coefficients[,4])
    }
    bt <- boot(data=data, statistic=bs_p, R=1000)

The class "numeric" values returned from each of these appears to be in exactly the same format, to my beginner's eye... but I'm guessing it isn't? I have also cleared the returned bt bootstrap object before running the next function, but that did not solve it. How could I best retrieve boot-strapped p-values? Thanks for any thoughts. (Running R 3.0.1 on Mac OSX.)

Answer 1

I am not sure if you can bootstrap p-values from lm model (but the solution is provided for that) . In your bs or bs_r function, you can remove d$ on the right hand side of fit since you already defined data d. Here is the example using mtcars data :

library(boot)

bs <- function(mtcars, i) {
  d <- mtcars[i,]
  fit <- lm (mpg~drat+wt, data=d)
  return(coef(fit))
}
bt <- boot(data=mtcars, statistic=bs, R=1000)

bt

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = mtcars, statistic = bs, R = 1000)


Bootstrap Statistics :
     original      bias    std. error
t1* 30.290370  0.54284222    7.494441
t2*  1.442491 -0.07260619    1.393801
t3* -4.782890 -0.09804271    1.000838

Here is the p-values for bootstrapped p-values from lm.

bs_r <- function(mtcars, i) {
  d <- mtcars[i,]
  fit <- lm (mpg~drat+wt, data=d)
  return(summary(fit)$coefficients[,4])
}
bt1 <- boot(data=mtcars, statistic=bs_r, R=1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = mtcars, statistic = bs_r, R = 1000)


Bootstrap Statistics :
        original       bias     std. error
t1* 2.737824e-04 4.020024e-03 0.0253248217
t2* 3.308544e-01 7.108738e-02 0.2960776146
t3* 1.589075e-06 5.405459e-05 0.0005540412