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