R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("bootstrap-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('bootstrap') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "Rainfall" > > ### * Rainfall > > flush(stderr()); flush(stdout()) > > ### Name: Rainfall > ### Title: Rainfall Data > ### Aliases: Rainfall > ### Keywords: datasets > > ### ** Examples > > str(Rainfall) Time-Series [1:106] from 1873 to 1978: 80 40 65 46 68 32 58 60 61 60 ... > plot(Rainfall) > > > > cleanEx(); ..nameEx <- "abcnon" > > ### * abcnon > > flush(stderr()); flush(stdout()) > > ### Name: abcnon > ### Title: Nonparametric ABC Confidence Limits > ### Aliases: abcnon > ### Keywords: nonparametric htest > > ### ** Examples > > # compute abc intervals for the mean > x <- rnorm(10) > theta <- function(p,x) {sum(p*x)/sum(p)} > results <- abcnon(x, theta) > # compute abc intervals for the correlation > x <- matrix(rnorm(20),ncol=2) > theta <- function(p, x) + { + x1m <- sum(p * x[, 1])/sum(p) + x2m <- sum(p * x[, 2])/sum(p) + num <- sum(p * (x[, 1] - x1m) * (x[, 2] - x2m)) + den <- sqrt(sum(p * (x[, 2] - x1m)^2) * + sum(p * (x[, 2] - x2m)^2)) + return(num/den) + } > results <- abcnon(x, theta) > > > > cleanEx(); ..nameEx <- "abcpar" > > ### * abcpar > > flush(stderr()); flush(stdout()) > > ### Name: abcpar > ### Title: Parametric ABC Confidence Limits > ### Aliases: abcpar > ### Keywords: htest > > ### ** Examples > > # binomial > # x is a p-vector of successes, n is a p-vector of > # number of trials > ## Not run: > ##D S <- matrix(0,nrow=p,ncol=p) > ##D S[row(S)==col(S)] <- x*(1-x/n) > ##D mu <- function(eta,n){n/(1+exp(eta))} > ##D etahat <- log(x/(n-x)) > ##D #suppose p=2 and we are interested in mu2-mu1 > ##D tt <- function(mu){mu[2]-mu[1]} > ##D x <- c(2,4); n <- c(12,12) > ##D a <- abcpar(x, tt, S, etahat,n) > ## End(Not run) > > > cleanEx(); ..nameEx <- "bcanon" > > ### * bcanon > > flush(stderr()); flush(stdout()) > > ### Name: bcanon > ### Title: Nonparametric BCa Confidence Limits > ### Aliases: bcanon > ### Keywords: nonparametric > > ### ** Examples > > # bca limits for the mean > # (this is for illustration; > # since "mean" is a built in function, > # bcanon(x,100,mean) would be simpler!) > x <- rnorm(20) > theta <- function(x){mean(x)} > results <- bcanon(x,100,theta) > > # To obtain bca limits for functions of more > # complex data structures, write theta > # so that its argument x is the set of observation > # numbers and simply pass as data to bcanon > # the vector 1,2,..n. > # For example, find bca limits for > # the correlation coefficient from a set of 15 data pairs: > xdata <- matrix(rnorm(30),ncol=2) > n <- 15 > theta <- function(x,xdata){ cor(xdata[x,1],xdata[x,2]) } > results <- bcanon(1:n,100,theta,xdata) > > > > cleanEx(); ..nameEx <- "bootpred" > > ### * bootpred > > flush(stderr()); flush(stdout()) > > ### Name: bootpred > ### Title: Bootstrap Estimates of Prediction Error > ### Aliases: bootpred > ### Keywords: nonparametric > > ### ** Examples > > # bootstrap prediction error estimation in least squares > # regression > x <- rnorm(85) > y <- 2*x +.5*rnorm(85) > theta.fit <- function(x,y){lsfit(x,y)} > theta.predict <- function(fit,x){ + cbind(1,x)%*%fit$coef + } > sq.err <- function(y,yhat) { (y-yhat)^2} > results <- bootpred(x,y,20,theta.fit,theta.predict, + err.meas=sq.err) > > # for a classification problem, a standard choice > # for err.meas would simply count up the > # classification errors: > miss.clas <- function(y,yhat){ 1*(yhat!=y)} > # with this specification, bootpred estimates > # misclassification rate > > > > cleanEx(); ..nameEx <- "bootstrap" > > ### * bootstrap > > flush(stderr()); flush(stdout()) > > ### Name: bootstrap > ### Title: Non-Parametric Bootstrapping > ### Aliases: bootstrap > ### Keywords: nonparametric > > ### ** Examples > > # 100 bootstraps of the sample mean > # (this is for illustration; since "mean" is a > # built in function, bootstrap(x,100,mean) would be simpler!) > x <- rnorm(20) > theta <- function(x){mean(x)} > > results <- bootstrap(x,100,theta) > > # as above, but also estimate the 95th percentile > # of the bootstrap dist'n of the mean, and > # its jackknife-after-bootstrap standard error > > perc95 <- function(x){quantile(x, .95)} > > > results <- bootstrap(x,100,theta, func=perc95) > > # To bootstrap functions of more complex data structures, > # write theta so that its argument x > # is the set of observation numbers > # and simply pass as data to bootstrap the vector 1,2,..n. > # For example, to bootstrap > # the correlation coefficient from a set of 15 data pairs: > xdata <- matrix(rnorm(30),ncol=2) > n <- 15 > theta <- function(x,xdata){ cor(xdata[x,1],xdata[x,2]) } > results <- bootstrap(1:n,20,theta,xdata) > > > > cleanEx(); ..nameEx <- "boott" > > ### * boott > > flush(stderr()); flush(stdout()) > > ### Name: boott > ### Title: Bootstrap-t Confidence Limits > ### Aliases: boott > ### Keywords: nonparametric > > ### ** Examples > > # estimated confidence points for the mean > x <- rchisq(20,1) > theta <- function(x){mean(x)} > results <- boott(x,theta) > # estimated confidence points for the mean, > # using variance-stabilization bootstrap-T method > results <- boott(x,theta,VS=TRUE) > results$confpoints # gives confidence points 0.001 0.01 0.025 0.05 0.1 0.5 0.9 [1,] 0.544524 0.5527472 0.6109023 0.6630976 0.7067375 1.083742 1.978057 0.95 0.975 0.99 0.999 [1,] 2.442318 3.030949 3.251848 3.315702 > # plot the estimated var stabilizing transformation > plot(results$theta,results$g) > # use standard formula for stand dev of mean > # rather than an inner bootstrap loop > sdmean <- function(x, ...) + {sqrt(var(x)/length(x))} > results <- boott(x,theta,sdfun=sdmean) > > # To bootstrap functions of more complex data structures, > # write theta so that its argument x > # is the set of observation numbers > # and simply pass as data to boot the vector 1,2,..n. > # For example, to bootstrap > # the correlation coefficient from a set of 15 data pairs: > xdata <- matrix(rnorm(30),ncol=2) > n <- 15 > theta <- function(x, xdata){ cor(xdata[x,1],xdata[x,2]) } > results <- boott(1:n,theta, xdata) > > > > cleanEx(); ..nameEx <- "cell" > > ### * cell > > flush(stderr()); flush(stdout()) > > ### Name: cell > ### Title: Cell Survival data > ### Aliases: cell > ### Keywords: datasets > > ### ** Examples > > plot(cell[,2:1],pch=c(rep(1,12),17,1), + col=c(rep("black",12),"red", "black"), + cex=c(rep(1,12), 2, 1)) > > > > cleanEx(); ..nameEx <- "cholost" > > ### * cholost > > flush(stderr()); flush(stdout()) > > ### Name: cholost > ### Title: The Cholostyramine Data > ### Aliases: cholost > ### Keywords: datasets > > ### ** Examples > > str(cholost) `data.frame': 164 obs. of 2 variables: $ z: num 0 27 71 95 0 28 71 95 0 29 ... $ y: num -5.25 -1.50 59.50 32.50 -7.25 ... > summary(cholost) z y Min. : 0.00 Min. :-23.000 1st Qu.: 27.00 1st Qu.: 7.125 Median : 69.50 Median : 30.500 Mean : 60.12 Mean : 32.806 3rd Qu.: 95.00 3rd Qu.: 54.562 Max. :100.00 Max. :113.250 > plot(y ~ z, data=cholost, xlab="Compliance", + ylab="Improvement") > abline(lm(y ~ z, data=cholost), col="red") > > > > cleanEx(); ..nameEx <- "crossval" > > ### * crossval > > flush(stderr()); flush(stdout()) > > ### Name: crossval > ### Title: K-fold Cross-Validation > ### Aliases: crossval > ### Keywords: nonparametric > > ### ** Examples > > # cross-validation of least squares regression > # note that crossval is not very efficient, and being a > # general purpose function, it does not use the > # Sherman-Morrison identity for this special case > x <- rnorm(85) > y <- 2*x +.5*rnorm(85) > theta.fit <- function(x,y){lsfit(x,y)} > theta.predict <- function(fit,x){ + cbind(1,x)%*%fit$coef + } > results <- crossval(x,y,theta.fit,theta.predict,ngroup=6) > > > > > cleanEx(); ..nameEx <- "diabetes" > > ### * diabetes > > flush(stderr()); flush(stdout()) > > ### Name: diabetes > ### Title: Blood Measurements on 43 Diabetic Children > ### Aliases: diabetes > ### Keywords: datasets > > ### ** Examples > > plot(logCpeptide ~ age, data=diabetes) > > > > cleanEx(); ..nameEx <- "hormone" > > ### * hormone > > flush(stderr()); flush(stdout()) > > ### Name: hormone > ### Title: Hormone Data from page 107 > ### Aliases: hormone > ### Keywords: datasets > > ### ** Examples > > str(hormone) `data.frame': 27 obs. of 3 variables: $ Lot : chr "A" "A" "A" "A" ... $ hrs : num 99 152 293 155 196 53 184 171 52 376 ... $ amount: num 25.8 20.5 14.3 23.2 20.6 31.1 20.9 20.9 30.4 16.3 ... > if(interactive())par(ask=TRUE) > with(hormone, stripchart(amount ~ Lot)) > with(hormone, plot(amount ~ hrs, pch=Lot)) > abline( lm(amount ~ hrs, data=hormone, col="red2")) Warning in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : extra arguments col are just disregarded. > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "jackknife" > > ### * jackknife > > flush(stderr()); flush(stdout()) > > ### Name: jackknife > ### Title: Jackknife Estimation > ### Aliases: jackknife > ### Keywords: nonparametric > > ### ** Examples > > # jackknife values for the sample mean > # (this is for illustration; # since "mean" is a > # built in function, jackknife(x,mean) would be simpler!) > x <- rnorm(20) > theta <- function(x){mean(x)} > > results <- jackknife(x,theta) > > # To jackknife functions of more complex data structures, > # write theta so that its argument x > # is the set of observation numbers > # and simply pass as data to jackknife the vector 1,2,..n. > # For example, to jackknife > # the correlation coefficient from a set of 15 data pairs: > > xdata <- matrix(rnorm(30),ncol=2) > n <- 15 > theta <- function(x,xdata){ cor(xdata[x,1],xdata[x,2]) } > results <- jackknife(1:n,theta,xdata) > > > > cleanEx(); ..nameEx <- "law" > > ### * law > > flush(stderr()); flush(stdout()) > > ### Name: law > ### Title: Law school data from Efron and Tibshirani > ### Aliases: law > ### Keywords: datasets > > ### ** Examples > > str(law) `data.frame': 15 obs. of 2 variables: $ LSAT: num 576 635 558 578 666 580 555 661 651 605 ... $ GPA : num 339 330 281 303 344 307 300 343 336 313 ... > if(interactive())par(ask=TRUE) > plot(law) > theta <- function(ind) cor(law[ind,1], law[ind,2]) > theta(1:15) # sample estimate [1] 0.7763745 > law.boot <- bootstrap(1:15, 2000, theta) > sd(law.boot$thetastar) # bootstrap standard error [1] 0.1318714 > hist(law.boot$thetastar) > # bootstrap t confidence limits for the correlation coefficient: > theta <- function(ind) cor(law[ind,1], law[ind,2]) > boott(1:15, theta, VS=FALSE)$confpoints 0.001 0.01 0.025 0.05 0.1 0.5 0.9 [1,] -3.769365 -0.8587222 -0.6335156 0.003080172 0.2288738 0.7702694 0.9529714 0.95 0.975 0.99 0.999 [1,] 1.004892 1.057899 1.176022 1.268214 > boott(1:15, theta, VS=TRUE)$confpoints 0.001 0.01 0.025 0.05 0.1 0.5 0.9 [1,] 0.1375004 0.2425394 0.3500606 0.398738 0.5040627 0.7574344 0.8953549 0.95 0.975 0.99 0.999 [1,] 0.9097642 0.9255485 0.9449852 0.9715406 > # Observe the difference! See page 162 of the book. > # abcnon(as.matrix(law), function(p,x) cov.wt(x, p, cor=TRUE)$cor[1,2] )$limits > # The above cannot be used, as the resampling vector can take negative values! > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "law82" > > ### * law82 > > flush(stderr()); flush(stdout()) > > ### Name: law82 > ### Title: Data for Universe of USA Law Schools > ### Aliases: law82 > ### Keywords: datasets > > ### ** Examples > > plot(law82[,2:3]) > cor(law82[,2:3]) LSAT GPA LSAT 1.0000000 0.7599979 GPA 0.7599979 1.0000000 > > > > cleanEx(); ..nameEx <- "lutenhorm" > > ### * lutenhorm > > flush(stderr()); flush(stdout()) > > ### Name: lutenhorm > ### Title: Luteinizing Hormone > ### Aliases: lutenhorm > ### Keywords: datasets > > ### ** Examples > > str(lutenhorm) `data.frame': 48 obs. of 5 variables: $ V1: num 1 2 3 4 5 6 7 8 9 10 ... $ V2: num 2.2 2.2 2.3 2 1.6 1.4 1.8 2.2 2.9 2.6 ... $ V3: num 5.5 4.5 5.1 5.5 5.7 5.1 4.3 4.8 5.6 5.9 ... $ V4: num 2.4 2.4 2.4 2.2 2.1 1.5 2.3 2.3 2.5 2 ... $ V5: num 4.3 4.6 4.7 4.1 4.1 5.2 5 4.4 4.2 5.1 ... > matplot(lutenhorm) > > > > cleanEx(); ..nameEx <- "mouse.c" > > ### * mouse.c > > flush(stderr()); flush(stdout()) > > ### Name: mouse.c > ### Title: Experiments with mouse > ### Aliases: mouse.c > ### Keywords: datasets > > ### ** Examples > > str(mouse.c) num [1:9] 52 104 146 10 50 31 40 27 46 > if(interactive())par(ask=TRUE) > stripchart(list(treatment=mouse.t, control=mouse.c)) > cat("bootstrapping the difference of means, treatment - control:\n") bootstrapping the difference of means, treatment - control: > cat("bootstrapping is done independently for the two groups\n") bootstrapping is done independently for the two groups > mouse.boot.c <- bootstrap(mouse.c, 2000, mean) > mouse.boot.t <- bootstrap(mouse.t, 2000, mean) > mouse.boot.diff <- mouse.boot.t$thetastar - mouse.boot.c$thetastar > hist(mouse.boot.diff) > abline(v=0, col="red2") > sd(mouse.boot.diff) [1] 26.97613 > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "mouse.t" > > ### * mouse.t > > flush(stderr()); flush(stdout()) > > ### Name: mouse.t > ### Title: Experiment with mouse > ### Aliases: mouse.t > ### Keywords: datasets > > ### ** Examples > > str(mouse.t) num [1:7] 94 197 16 38 99 141 23 > stripchart(list(treatment=mouse.t, control=mouse.c)) > > > > cleanEx(); ..nameEx <- "patch" > > ### * patch > > flush(stderr()); flush(stdout()) > > ### Name: patch > ### Title: The Patch Data > ### Aliases: patch > ### Keywords: datasets > > ### ** Examples > > str(patch) `data.frame': 8 obs. of 6 variables: $ subject : int 1 2 3 4 5 6 7 8 $ placebo : num 9243 9671 11792 13357 9055 ... $ oldpatch: num 17649 12013 19979 21816 13850 ... $ newpatch: num 16449 14614 17274 23798 12560 ... $ z : num 8406 2342 8187 8459 4795 ... $ y : num -1200 2601 -2705 1982 -1290 ... > theta <- function(ind){ + Y <- patch[ind,"y"] + Z <- patch[ind,"z"] + mean(Y)/mean(Z) } > patch.boot <- bootstrap(1:8, 2000, theta) > names(patch.boot) [1] "thetastar" "func.thetastar" "jack.boot.val" "jack.boot.se" [5] "call" > hist(patch.boot$thetastar) > abline(v=c(-0.2, 0.2), col="red2") > theta(1:8) #sample plug-in estimator [1] -0.0713061 > abline(v=theta(1:8) , col="blue") > # The bootstrap bias estimate: > mean(patch.boot$thetastar) - theta(1:8) [1] 0.003553931 > sd(patch.boot$thetastar) # bootstrapped standard error [1] 0.09971667 > > > > cleanEx(); ..nameEx <- "scor" > > ### * scor > > flush(stderr()); flush(stdout()) > > ### Name: scor > ### Title: Open/Closed Book Examination Data > ### Aliases: scor > ### Keywords: datasets > > ### ** Examples > > str(scor) `data.frame': 88 obs. of 5 variables: $ mec: num 77 63 75 55 63 53 51 59 62 64 ... $ vec: num 82 78 73 72 63 61 67 70 60 72 ... $ alg: num 67 80 71 63 65 72 65 68 58 60 ... $ ana: num 67 70 66 70 70 64 65 62 62 62 ... $ sta: num 81 81 81 68 63 73 68 56 70 45 ... > if(interactive())par(ask=TRUE) > plot(scor) > # The parameter of interest (theta) is the fraction of variance explained > # by the first principal component. > # For principal components analysis svd is better numerically than > # eigen-decomposistion, but for bootstrapping the later is MUCH faster. > theta <- function(ind) { + vals <- eigen(var(scor[ind,]), symmetric=TRUE, only.values=TRUE)$values + vals[1] / sum(vals) } > scor.boot <- bootstrap(1:88, 500, theta) > sd(scor.boot$thetastar) # bootstrap standard error [1] 0.04948759 > hist(scor.boot$thetastar) > abline(v=theta(1:88), col="red2") > abline(v=mean(scor.boot$thetastar), col="blue") > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "spatial" > > ### * spatial > > flush(stderr()); flush(stdout()) > > ### Name: spatial > ### Title: Spatial Test Data > ### Aliases: spatial > ### Keywords: datasets > > ### ** Examples > > str(spatial) `data.frame': 26 obs. of 2 variables: $ A: num 48 36 20 29 42 42 20 42 22 41 ... $ B: num 42 33 16 39 38 36 15 33 20 43 ... > plot(spatial) > abline(0,1, col="red2") > > > > cleanEx(); ..nameEx <- "stamp" > > ### * stamp > > flush(stderr()); flush(stdout()) > > ### Name: stamp > ### Title: Data on Thickness of Stamps > ### Aliases: stamp > ### Keywords: datasets > > ### ** Examples > > summary(stamp) Thickness Min. :0.06000 1st Qu.:0.07500 Median :0.08000 Mean :0.08602 3rd Qu.:0.09800 Max. :0.13100 > with(stamp, {hist(Thickness); + plot(density(Thickness), add=TRUE)}) Warning: parameter "add" could not be set in high-level plot() function Warning: parameter "add" could not be set in high-level plot() function Warning: parameter "add" could not be set in high-level plot() function Warning: parameter "add" could not be set in high-level plot() function Warning: parameter "add" could not be set in high-level plot() function Warning: parameter "add" could not be set in high-level plot() function > > > > cleanEx(); ..nameEx <- "tooth" > > ### * tooth > > flush(stderr()); flush(stdout()) > > ### Name: tooth > ### Title: Tooth Strength Data > ### Aliases: tooth > ### Keywords: datasets > > ### ** Examples > > str(tooth) `data.frame': 13 obs. of 6 variables: $ patient : int 1 2 3 4 5 6 7 8 9 10 ... $ D1 : num -5.29 -5.94 -5.61 -5.41 -5.20 ... $ D2 : num 10.09 10.00 10.18 10.13 8.84 ... $ E1 : num 12.3 11.4 11.8 12.1 10.7 ... $ E2 : num 13.1 13.0 13.2 12.8 11.7 ... $ strength: num 36.0 35.5 35.4 36.0 34.6 ... > mod.easy <- lm(strength ~ E1+E2, data=tooth) > mod.diffi <- lm(strength ~ D1+D2, data=tooth) > summary(mod.easy) Call: lm(formula = strength ~ E1 + E2, data = tooth) Residuals: Min 1Q Median 3Q Max -1.1681 -0.1783 0.1203 0.4313 0.5846 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 27.66657 5.02179 5.509 0.000258 *** E1 0.59761 0.39095 1.529 0.157355 E2 0.03716 0.55306 0.067 0.947754 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5593 on 10 degrees of freedom Multiple R-Squared: 0.3336, Adjusted R-squared: 0.2003 F-statistic: 2.503 on 2 and 10 DF, p-value: 0.1314 > summary(mod.diffi) Call: lm(formula = strength ~ D1 + D2, data = tooth) Residuals: Min 1Q Median 3Q Max -1.0843 -0.2951 0.1411 0.4049 0.5396 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 30.4915 3.8726 7.874 1.35e-05 *** D1 0.6991 0.4337 1.612 0.1380 D2 0.8637 0.3597 2.401 0.0373 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.5254 on 10 degrees of freedom Multiple R-Squared: 0.412, Adjusted R-squared: 0.2944 F-statistic: 3.504 on 2 and 10 DF, p-value: 0.07028 > if(interactive())par(ask=TRUE) > theta <- function(ind) { + easy <- lm(strength ~ E1+E2, data=tooth, subset=ind) + diffi<- lm(strength ~ D1+D2, data=tooth, subset=ind) + (sum(resid(easy)^2) - sum(resid(diffi)^2))/13 } > tooth.boot <- bootstrap(1:13, 2000, theta) > hist(tooth.boot$thetastar) > abline(v=0, col="red2") > qqnorm(tooth.boot$thetastar) > qqline(tooth.boot$thetastar, col="red2") > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > ### *