Tuesday, June 11, 2013

Sample R code for bootstrapping residues for regression

y = rnorm(100) + rep(c(0, 1), each = 50)
x = rep(c(1,2), each = 50)

my = mean(y)
res = lm(y~x)$resid

# Log LR test
lr.obs = sum((y-my)^2) - sum(res^2) # Can replace with "F" like statistic

lr.star = c()
for(boot in 1:1000){
y.star = my + sample(res, replace = T)
res.star = lm(y.star~x)$resid
lr.star[boot] = sum((y.star-my)^2) - sum(res.star^2)
}

pval = mean(lr.star>= lr.obs)



This sample code was provided by Dr. Dipen Sangurdekar.  Boostrapping residue of the empirical data is considered to be more reasonable estimation of p-values.

No comments:

Post a Comment